# Multi-Spectral Reflection Matrix for Ultra-Fast 3D Label-Free Microscopy

Paul Balondrade, Victor Barolle, Nicolas Guigui, Emeric Auriant, Nathan

Rougier, Claude Boccara, Mathias Fink, Alexandre Aubry

► To cite this version:

Paul Balondrade, Victor Barolle, Nicolas Guigui, Emeric Auriant, Nathan Rougier, et al. Multi-Spectral Reflection Matrix for Ultra-Fast 3D Label-Free Microscopy. Nature Photonics, 2024, 10.1038/s41566-024-01479-y. hal-04211838v3

HAL Id: hal-04211838

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

Submitted on 11 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Multi-Spectral Reflection Matrix for

Ultra-Fast 3D Label-Free Microscopy

Paul Balondrade†, Victor Barolle†, Nicolas Guigui, Emeric Auriant,

Nathan Rougier, Claude Boccara, Mathias Fink, and Alexandre Aubry*

Institut Langevin, ESPCI Paris,

PSL University, CNRS, 75005 Paris, France

† These authors equally contributed to this work.

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

(Dated: September 11, 2024)
Abstract

Label-free microscopy exploits light scattering to obtain a three-
dimensional image of biological tissues. However, light propagation is af-
fected by aberrations and multiple scattering, which drastically degrade the
image quality and limit the penetration depth. Multi-conjugate adaptive op-
tics and time-gated matrix approaches have been developed to compensate
for aberrations but the associated frame rate is extremely limited for 3D
imaging. Here we develop a multi-spectral matrix approach to solve these
fundamental problems. Based on a sparse illumination scheme and an inter-
ferometric measurement of the reflected wave-field at multiple wavelengths,
the focusing process can be optimized in post-processing for any voxel by ad-
dressing independently each frequency component of the reflection matrix.
A proof-of-concept experiment demonstrates the three-dimensional image of
an opaque human cornea over a 0.1 mm³-field-of-view at a 290 nm-resolution
and a 1 Hz-frame rate. This work paves the way towards a fully-digital
microscope allowing real-time, in-vivo, quantitative and deep inspection of
tissues.
## 10 Introduction

Imaging of thick scattering tissues remains the greatest challenge in label-free microscopy$^{1-3}$. On the one hand, short-scale inhomogeneities of the refractive index backscatter light and the reflected wave-field can be leveraged to provide a structural image of the sample. On the other hand, larger-scale inhomogeneities give rise to forward multiple scattering events that distort the incident and reflected wave-fronts. This phenomenon, known as aberrations, leads to a drastic degradation of resolution and contrast at depths greater than the scattering mean free path $\ell_s$ ($\sim 100$ µm in biological tissues).

To circumvent this issue, adaptive optics (AO) has been transposed from astronomy to microscopy for the last twenty years$^4$. The basic idea is to compensate for wave distortions either by a direct sampling of the wave-field generated by a guide star or by an indirect metric optimization of the image. Unfortunately, AO correction is limited to a finite area, the so-called isoplanatic patch, the area over which aberrations can be considered spatially invariant. This problem becomes particularly important for deep imaging, where each isoplanatic patch reduces to a speckle grain at depths larger than the transport mean free path $\ell_t$ ($\sim 1$ mm in biological tissues). Multi-conjugate AO could increase the corrected field-of-view$^5$, but this would be at the price of a much more complex optical setup and an extremely long optimization process$^6$.

More recently, following seminal works that proposed post-processing computational strategies for AO$^{7-11}$, a reflection matrix approach has been developed for deep imaging$^{12-18}$. The basic idea is to illuminate the sample by a set of input wave-fronts and record via interferometry the reflected wave-front on a camera. Once this reflection matrix is measured, a set of matrix operations can be applied in order to perform a local compensation of aberrations and restore a diffraction-limited resolution for each pixel of the field-of-view. Nevertheless, the existing approaches suffer from several limitations. In most experimental works$^{12,13,15-18}$,
38 the reflection matrix is time-gated around the ballistic time as usually performed
39 in time-domain OCT¹⁹. Such a measurement has one main advantage since it en-
40 ables the temporal filtering of most multiply-scattered photons²⁰. However, it also
41 suffers from two strong drawbacks. First, time-gating means that a large part of
42 the information on the medium is discarded: Only the weakly distorted paths are
43 recorded and can be compensated by a spatial phase modulation of the incident
44 and reflected wave-fronts. Second, volumetric imaging can only be obtained by a
45 mechanical axial scanning of the sample, which limits the frame rate *Fps* to, at
46 best, 10⁶ voxels.s⁻¹ for a high quality correction over millimetric FOVs.

47 To go beyond, an acquisition of a spectral reflection matrix is required in or-
48 der to capture all the information required for the three-dimensional imaging of
49 a sample. In recent works²¹,²², the spatio-temporal degrees of freedom exhibited
50 by the reflection matrix have been exploited for tailoring dispersive focusing laws.
51 However, the acquisition rate was slow (Fps ~ 10³ voxels.s⁻¹) because the num-
52 ber of input wave-fronts scaled as the number of voxels in the image. Moreover,
53 the experimental demonstration was limited to the imaging of a resolution target
54 through a scattering medium²¹,²² or a sparse medium made of colloidal particles²².
55 In this paper, we go beyond an academic proof-of-concept and address the ex-
56 tremely challenging case of ultra-fast 3D imaging of biological tissues themselves
57 (nerves, cells, collagen, extracellular matrix etc.). In particular, we will show how
58 the number of input wave-fronts can be drastically decreased by deterministic fo-
59 cusing operations applied to the reflection matrix guided by a self-portrait of the
60 focusing process.

61 To that aim, we report on a measurement of the multi-spectral reflection ma-
62 trix at a much higher frame rate (Fps ~ 5.10⁹ voxels.s⁻¹), with a 3D imaging
63 demonstration on an *ex-vivo* opaque cornea at a resolution of 0.29 µm and 0.5
64 µm in the transverse and axial directions, respectively. The experimental set up
65 combines a Fourier-domain full-field OCT (FD-FF-OCT) setup²³⁻²⁵ with a coher-
ent multi-illumination scheme. Capable of recording a polychromatic reflection
matrix of $10^{10}$ coefficients in less than 1 s with an ultra-fast camera, this device is
fully compatible with future *in-vivo* applications. As in FD-FF-OCT, a frequency
Fourier transform and numerical refocusing can provide a 3D image of the sample
for each incident wave-front<sup>23-25</sup> but, as expected, multiple scattering is shown to
strongly hamper the imaging process. A coherent compound of images obtained
for each illumination in post-processing can then provide a digital confocal image
but its resolution and contrast are drastically affected by sample-induced aberra-
tions. Interestingly, reflection matrix imaging (RMI) can go beyond by decoupling
input and output focusing points at each time-of-flight. A focused and time-gated
reflection matrix is synthesized. It contains the impulse responses between virtual
sources and detectors conjugated with each voxel inside the sample. While pre-
vious works only considered focusing points at the same depth<sup>13-16</sup>, we show here
that their axial scan gives access to a self-portrait of light propagation inside the
sample. The focusing process can then be guided by this reflection point spread
function, thereby allowing an optimized compensation of defocus aberrations at
each depth of the sample, even in presence of multiple scattering. Finally, trans-
verse aberrations can be tackled by means of a local analysis of wave distortions,
as initially proposed with ultrasound<sup>26</sup> and recently transposed to optics<sup>18</sup>. How-
ever, while those previous works implied a complex multi-scale analysis of wave
distortions, the prior compensation of defocus aberrations here enables a direct
and local correction of transverse aberrations. A digital clearing of refractive in-
dex heterogeneities is thus applied and a three-dimensional image of the sample
is obtained with an optimized contrast and close-to-ideal resolution throughout
the volume. The ability of RMI in addressing forward multiple scattering paths
is finally investigated by considering the imaging of a resolution target through a
more opaque region of the cornea.
<figure><img src="image_2.png" /><figcaption>FIG. 1. **Measuring the multi-spectral reflection matrix.** **a,** A wavelength swept light source illuminates a Linnik interferometer through a collimator, two scanning mirrors and a lens (L1) that allows a raster scanning of the focal spot in the MO pupil planes ($\mathbf{u}_{\text{in}}$) in each arm ($\mathbf{a}_1$). The sample placed in the focal plane of the first MO (MO₁, NA=0.8) is thus illuminated by a set of plane waves at each frequency of the light source bandwidth ($\mathbf{a}_2$). The backscattered wave field is collected through the same MO, focused by means of a second lens L2 on the surface of a CMOS camera where it interferes with a reference beam ($\mathbf{a}_3$). The latter beam results from the reflection of the same incident wave-fronts by a reference mirror placed in the focal plane of the second MO (MO₂, NA=0.8). **b,** At each frequency $\omega$, for each input wave-front $\mathbf{u}_{\text{in}}$, the interferogram $I(\mathbf{s}_{\text{out}})$ ($\mathbf{b}_1$) recorded by each pixel $\mathbf{s}_{\text{out}}$ of the camera provides one column of the spectral reflection matrix $\mathbf{R}_{\text{su}}(\omega) = [R(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega)]$ ($\mathbf{b}_2$).</figcaption></figure>

93 **Results**

94 **Recording the Multi-Spectral Reflection Matrix.**

95 3D matrix imaging is based on the measurement of a multi-spectral reflection 96 matrix from the scattering sample. The experimental setup and procedure are 97 described in Fig. 1 (see Methods and Supplementary Figure S1). Inspired by FD-98 FF-OCT²⁷, it simply consists in a Linnik interferometer (Fig. 1a). In the first arm, 99 a reference mirror is placed in the focal plane of a microscope objective (MO). The
100 second arm contains the scattering sample to be imaged through an identical MO.
101 This interferometer is illuminated by a swept source through two scanning mirrors
102 and a lens that allows a raster scanning of the focused beam in the MO pupil
103 planes (Fig. 1a₁). The sample and reference mirror are thus illuminated by a set
104 of plane waves at each frequency of the light source bandwidth (Fig. 1a₂). The
105 reflected waves are collected through the same MOs and, ultimately, interfere on a
106 camera conjugated with the focal plane. For each input wave-front of frequency ω
107 and coordinate **u**<sub>in</sub> in the pupil plane, an interferogram *I*(*s*<sub>out</sub>, **u**<sub>in</sub>, ω) is recorded by
108 the camera whose sensors are identified by their transverse position **s**<sub>out</sub> (Fig. 1a₃).
109 Each interferogram provides one column of the multi-spectral reflection matrix
110 **R**<sub>su</sub>(ω) = [*R*(*s*<sub>out</sub>, **u**<sub>in</sub>, ω)] (see Methods and Fig. 1b). The subscripts associated
111 with **R** here denote the output and input basis in which the reflection matrix
112 is defined. Indeed, this matrix will be projected between different bases in this
113 work: (i) the camera sensor basis (s); (ii) the plane wave basis (u) and; (iii) the
114 focused basis (ρ), in which the image is built. The coefficients of the multi-spectral
115 reflection matrix, *R*(*s*<sub>out</sub>, **u**<sub>in</sub>, ω), are noted as functions of the coordinates in the
116 planes associated with each basis. The same convention will be considered for the
117 other operators considered in this paper.

118 In the opaque cornea experiment, the reflection matrix **R** is measured with
119 $N_{\text{in}} = 177$ plane waves, corresponding to a full scan of the immersion MO pupil
120 (NA=0.8, refractive index $n_0 = 1.33$). The interferograms are recorded by $N_{\text{out}} =$
121 $1024^2$ pixels of the camera, corresponding to an output FOV of $\Omega_{\text{out}} \times \Omega_{\text{out}} = 297 \times$
122 $297 \ \mu\text{m}^2$, with a spatial sampling $\delta\rho_{\text{out}} = 290 \ \text{nm}$. Finally, $N_{\omega} = 201$ independent
123 frequencies are used to probe the sample within the frequency bandwidth [$800,875$] 
124 nm of the light source. All the information about the sample is thus contained
125 in the $2.10^9$ coefficients acquired in 1.4 s. In the following, we show how to post-
126 process this wealth of optical data to build a 3D highly-contrasted image of the
127 cornea at a diffraction-limited resolution.
128 Ultra-fast Three-Dimensional Imaging.

To that aim, the most direct path is to perform, a Fourier transform over frequency $\omega$ of the back-scattered wave-field recorded for one plane-wave illumination$^{23}$: This is the principle of full-field swept-source OCT which provides an image whose axial dimension is dictated by photons' times-of-flight (Supplementary Figure S2). The resulting image is, however, completely blurred without any connection with the sample reflectivity. Indeed, a high NA implies a very restricted depth-of-field which is prohibitory for 3D imaging$^{28}$: $\delta z_f \sim 2n_0\lambda/NA^2 \sim 3.5$ µm, with $\lambda$ the central wavelength. A prior numerical focusing of the wave-field recorded by the camera shall be performed at each depth $z$ of the sample. This is the principle of the holoscope developed by Hillmann et al. about a decade ago$^{10}$.

This numerical focusing process is performed by means of Fresnel propagators. For this purpose, the multi-spectral reflection matrix is projected at output in the focused basis:

$$ \mathbf{R}_{\rho u}(z, \omega) = \mathbf{F}_{\rho s}^*(z, \omega) \times \mathbf{R}_{s u}(\omega) \quad (1) $$

where the symbols $*$ and $\times$ stand for the phase conjugation operation and the standard matrix product, respectively. $\mathbf{F}_{\rho s}(z, \omega) = [F(\rho_{\text{out}}, \mathbf{s}_{\text{out}}, z, \omega)]$ is the Fresnel operator that describes free-space propagation from the camera ($\mathbf{s}$) to any focal plane ($\rho$) located at expected depth $z$ in the sample (Methods, Eq. 22). Each frequency component of $\mathbf{R}_{\rho u}(z, \omega)$ should then be recombined in order to time gate the singly-scattered photons. In practice, an inverse Fourier transform over frequency $\omega$ is performed and yields a $\mathbf{R}$-matrix as a function of photon's time-of-flight $t$:

$$ \mathbf{R}_{\rho u}(z, t) = \int d\omega \mathbf{R}_{\rho u}(z, \omega) e^{j\omega t}. \quad (2) $$

The coherence volume is defined as the ensemble of points that contribute to the singly-scattered wave-field at a given time-of-flight $t$. In an homogeneous medium of refractive index $n$, this volume is, in first approximation, a thin slice located
153 at a depth $z_t = c_0t/(2n)$, with $c_0$ the light velocity in vacuum. Its thickness is 154 governed by the light source bandwidth $\Delta\lambda$: $\delta z_t \sim \lambda^2/(2n_0\Delta\lambda) \sim 3.5$ µm. When 155 the focusing plane and the coherence volume coincide (Fig. 2a₁), an holoscopic 156 image of the sample, $I_H$, can be obtained for each input wave-front $\mathbf{u}_{\text{in}}$ (Fig. 2a₁):

$$I_H(\mathbf{r}_t, \mathbf{u}_{\text{in}}) = R(\rho_{\text{out}}, \mathbf{u}_{\text{in}}, z_t, t). \qquad (3)$$

157 with $\mathbf{r}_t = (\rho_{\text{out}}, z_t)$. In practice, an exact matching between the focusing plane 158 and coherence volume is difficult to obtain especially for deep imaging (i.e low 159 single-to-multiple scattering ratio). We will describe further how matrix imaging 160 can provide a robust observable for this fine tuning.

161 Figures 2b₁, c₁ and d₁ display longitudinal and transverse cross-sections of the 162 cornea obtained for a normal incident plane wave (see also Supplementary Movies 163 1 and 2). Although this holoscopic image can be obtained at a very high frame 164 rate (Fps ~ $10^{11}$ voxels/s), it also exhibits a speckle-like feature. Indeed, multi- 165 ple scattering events taking place ahead of the coherence volume at each time t 166 can pollute the image. Multiply-scattered waves generate a random speckle noise 167 without any connection with the medium reflectivity. To remove it, a naive strat- 168 egy is to sum the intensity of the holoscopic images obtained for each illumination 169 $\mathbf{u}_{\text{in}}$. Such an incoherent compound tends to smooth out the speckle noise but 170 the resulting image still exhibits an extremely low contrast due to the multiple 171 scattering background (see Supplementary Fig. S2). To get rid of it, the single-to- 172 multiple scattering ratio shall be increased²⁰. For this purpose, a spatial filtering 173 of multiply-scattered photons can be performed by means of a confocal filter. Nev- 174 ertheless, this operation is extremely sensitive to the focusing quality inside the 175 sample. A prior optimization of the focusing process is thus needed.

176 **Digital confocal imaging.**

177 To that aim, the dual reflection matrix is projected in the focused basis at input.
<figure><img src="image_3.png" /><figcaption>FIG. 2. **From Holoscopy to Matrix Imaging.** **a,** Imaging Methods. **a**₁, Holoscopy: The sample is illuminated by a plane wave (in green) and an image is produced by spatio-temporal focusing of the back-scattered wave-field on each voxel (ρ out , z) mapping the sample (red). **a**₂, Digital confocal microscopy (DCM): The sample is illuminated by a set of plane waves (in green) and a focused refection matrix **R** ρρ (z t ) is built by numerical focusing. A 3D confocal image is deduced from the diagonal elements of **R** ρρ (z t ) at each depth z t . **a**₃, Reflection matrix imaging (RMI): A local compensation of wave distortions is performed for each voxel. **b,** B-scan image showing one longitudinal section of the cornea reflectivity. **c-d,** En-face image of the cornea at z t = 105 µm and 230 µm, respectively [scale bar: 75 µm]. In panels (a)-(d), subscripts 1, 2 and 3 stand for holoscopy, DCM and RMI, respectively. In panel c₃, blue and red arrows design some corneal nerves and keratocytes, respectively; in panel d₃, stroma striae are highlighted by green dashed lines. **e,** RPSF at depth z t = 105 µm [scale bar: 3 µm] for DCM (e₁-f₁) and RMI (e₂-f₂) images. The color bar is in a square root scale. The radial evolutions of these RPSFs are compared in panel (d) [DCM: blue; RMI: orange]. **f,** Same as in panel (e) but at depth z t = 230 µm.</figcaption></figure>
179 Mathematically, it simply consists in a numerical input focusing of $\mathbf{R}_{\rho u}$ using the 180 Fresnel propagator $\mathbf{F}_{\rho u}$ that describes free space propagation from the MO pupil 181 plane and the focal plane at depth $z$ (Methods, Eq. 23):

$$ \mathbf{R}_{\rho\rho}(z, t) = \int d\omega \mathbf{R}_{\rho u}(z, \omega) \times \mathbf{F}_{\rho u}^{\dagger}(z, \omega) e^{-j\omega t}. \quad (4) $$

182 where the symbol $\dagger$ stands for transpose conjugate. Expressed in the focused 183 basis, the reflection matrix $\mathbf{R}_{\rho\rho}(z, t)$ contains the responses at each time-of-flight 184 $t$ between virtual sensors of expected positions $(\rho_{in}, z)$ and $(\rho_{out}, z)$.

185 The focused $\mathbf{R}$-matrix is equivalent to the time-gated reflection matrix consid- 186 ered in previous studies for RMI$^{12–18}$ except that we now have at our disposal a 187 supplementary degree of freedom: the parameter $z$ that controls the axial position 188 of the focusing plane. A raw confocal image $\mathbf{I}_C$ can be built by considering the 189 diagonal elements of $\mathbf{R}_{\rho\rho}$ ($\rho_{in} = \rho_{out}$):

$$ I_C(\mathbf{r}_t, z) = R(\rho, \rho, z, t). \quad (5) $$

190 Figure 3c shows the en-face image obtained at a given time-of-flight $t$ for different 191 values of $z$. Qualitatively, we see that the image quality strongly depends on the 192 relative position between the coherence volume and the focusing plane. Here the 193 presence of a highly reflecting structure, a corneal nerve, allows us to determine the 194 parameter $z$ that allows to match the focusing plane with the coherence volume.

195 **Self-portrait of the focusing process.**

196 A more quantitative and robust observable is provided by the off-diagonal co- 197 efficients of $\mathbf{R}_{\rho\rho}(z, t)$ that enable to probe the focusing quality at any voxel. More 198 precisely, this can be done by investigating the reflection point spread function
¹⁹⁹ (RPSF) defined as follows:

$$RPSF(\Delta\rho, \rho, z, t) = |R(\rho - \Delta\rho/2, \rho + \Delta\rho/2, z, t)|^2, \quad (6)$$

200 This quantity derived from the off-diagonal coefficients of $R_{\rho\rho}$, quantifies the focusing quality for each point $r_t = (\rho, z_t)$. For a medium of random reflectivity 201 and under a local isoplanatic assumption, its ensemble average actually scales as²⁹ 202 (Supplementary Section S5):

$$\langle RPSF(\Delta\rho, \rho, z, t) \rangle \propto |h_{\text{in}}|^2 \stackrel{\Delta\rho}{\otimes} |h_{\text{out}}|^2 (\Delta\rho, \rho, z, t) \quad (7)$$

204 where the symbols $\langle \cdots \rangle$ and $\otimes$ stand for ensemble average and convolution product, respectively. $h_{\text{in/out}}(\Delta\rho, \rho, z, t)$ is the spatial distribution of the input/output 205 PSF along the de-scanned coordinate $\Delta\rho$ in the coherence plane at $z_t$ when trying 206 to focus at point $(\rho, z)$. 208

209 The RPSF can thus provide a self-portrait of the focusing process inside the 210 cornea. Figure 3a shows the evolution of the laterally-averaged RPSF for a given 211 time $t$ as a function of the parameter $z$ in the Fresnel propagator. As expected, 212 the focusing plane and coherence volume coincide when the full width at half 213 maximum (FWHM) $\delta\rho$ of the RPSF is minimized (Fig. 3b), i.e for a defocus 214 distance $\Delta z = z - z_t = 0$ (Fig. 3a₂). The estimated defocus is roughly constant 215 over the whole thickness of the cornea. This proves that the effective index of the 216 cornea is actually very close to the refractive index $n_0$ used in our propagation 217 model (see Supplementary Section S4). 218

218 Figures 2b₂, c₂ and d₂ displays longitudinal and transverse cross-sections of the 219 confocal image obtained after tuning the coherence volume and focusing plane at 220 any depth (see also Supplementary Movies 1 and 2). The resolution and contrast 221 are much better than the incoherent compound image (Supplementary Fig. S2). In 222 particular, the axial resolution $\delta z_c$ of the digital confocal image is dictated by the
<figure><img src="image_4.png" /><figcaption>FIG. 3. **Auto-focusing process guided by the reflection point spread function.**
**a** Evolution of the RPSF versus the defocus distance $\Delta z$ for a fixed coherence volume. Its transverse distribution is shown for several values of defocus ($\mathbf{a}_1$)-(a3). The evolution of its radial average is displayed in panel (a4). The color bar in a square root scale. **b.** Relative position of the focusing plane (dash-dotted line) and coherence volume (red layer) for the different values of defocus $\Delta z$ displayed in (a). **c.** En-face confocal image and zoom on a nerve. In each panel, the subscripts 1, 2 and 3 stand for defocus distances $\Delta z = -4$, 0, and +4 $\mu$m. The considered coherence volume is located at the effective depth $z_t = 140 \mu$m in the cornea.</figcaption></figure>

223 high numerical aperture of the microscope objective²⁸: $\delta z_c \sim n_0 \lambda / (4NA^2) \sim 0.5$

224 $\mu$m. However, the image quality remains perfectible. Indeed, the RPSF still

225 spreads well beyond the theoretical resolution cell with $\delta\rho \sim 2.3 \mu$m at $z_t = 230$

226 $\mu$m (Fig. 2f₁). The Strehl ratio associated with this RPSF is extremely low ($S \sim 4.5 \times 10^{-3}$, Methods, Eq. 25) and the image contrast is relatively weak ($\chi \sim 4$ dB,
228 Methods, Eq. 26). These imperfections are the result of transverse aberrations
229 originating from the lateral fluctuations of the optical index $n(\rho, z)$ in the cornea.
230 To demonstrate this last assertion, the transverse evolution of the focusing process
231 can be investigated by a local assessment of the focusing quality (Methods, Eq. 24).
232 A map of local RPSFs is displayed in Fig. 4e₁. Although a prior compensation of
233 defocus aberrations provides a correct focusing quality over the whole thickness
234 of the cornea on average, the local RPSFs exhibit important fluctuations across-
235 the field-of-view. This observation is a manifestation of the 3D distribution of
236 the optical index $n(\mathbf{r})$ inside the cornea. This anisoplanic feature requires a local
237 compensation of aberrations as we will see below.

238 Local Compensation of Wave Distortions.

239 By considering the focused reflection matrix at each depth $z_t$ (Fig. 4a), a lo-
240 cal compensation of transverse aberrations can be performed using the distortion
241 matrix concept¹⁵. To that aim, the distorted component of the reflected/incident
242 wave-front is extracted for each input/output focusing point $\rho_{in/out}$ (Methods,
243 Eq. 28). Wave distortions exhibit local correlations that are a manifestation of the
244 shift-shift memory effect³⁰ characteristic of anisotropic scattering in the cornea
245 (Fig. 4b). Over an isoplanatic patch, each distorted wave-field corresponds to the
246 diffraction pattern of each virtual source ($\rho_{in/out}$) modulated by the local aberra-
247 tion transmittance that accounts for the long-scale heterogeneities of the refractive
248 index between the cornea surface and the coherence plane (Supplementary Section
249 S6). Each virtual source is spatially incoherent due to the random reflectivity of
250 the medium. The idea is to smartly combine each of them to synthesize a co-
251 herent guide star and estimate an aberration phase law independently from the
252 sample reflectivity. In practice, this is done by applying an iterative phase re-
253 versal (IPR) process²⁶ around each voxel ($\rho, z_t$) (Methods, Eq. 35). The result is
254 an estimation of the local aberration transmittance in the pupil plane¹⁸ (Fig. 4c).
255 The same operation can be repeated for each voxel of the medium. The set of
<figure><img src="image_5.png" /><figcaption>FIG. 4. **Matrix Compensation of Transverse Aberrations.** a, The digitally-refocused **R**-matrix contains the set of impulse responses $R(\rho_{\text{in}}, \rho_{\text{out}}, z_t)$ between an array of point sources $\rho_{\text{in}}$ and detectors $\rho_{\text{out}}$ lying in the coherence plane at depth $z_t$. b, Wave distortions are then isolated by projecting the optical data in the input or output pupil plane. Seen from the focal plane, this operation amounts to a de-scan of each focal spot. c, An iterative phase reversal process applied to those wave distortions provides an estimation of local aberration transmittances by recombining each focal spot into a virtual guide star over reduced spatial windows. d, The resulting maps of input ($d_1$) and output ($d_2$) aberration phase laws, $\phi_{\text{in/out}}$, are shown for depth $z_t = 140$ µm. e, The corresponding map of RPSFs is shown in amplitude before ($e_1$) and after ($e_2$) the compensation of transverse aberrations at the same depth. Each RPSF is displayed over a de-scan area of 6 × 6 µm$^2$. The color bar is in a square root scale.</figcaption></figure>

<footnote n="256">256 extracted aberration phase laws form the input and output aberration phase matrices, $\mathbf{\Phi}_{\text{in/out}}(z_t) = [\phi_{\text{in/out}}(\mathbf{u}_{\text{in/out}}, \rho, z_t)]$, between the pupil plane ($\mathbf{u}_{\text{in/out}}$) and the medium voxels $(\rho, z_t)^{26}$.</footnote>

<footnote n="259">259 A crucial parameter is the spatial extent $L$ of the spatial window over which</footnote>
260 correlations between wave distortions shall be investigated (Methods, Eq. 33).
261 This choice actually governs the transverse resolution of $\mathbf{\Phi}_{\text{in/out}}(z_t)$ and is subject
to the following dilemma: On the one hand, the spatial window should be as small
262 as possible to grasp the rapid variations of the aberration phase laws across the
263 field of view; on the other hand, these areas should be large enough to encompass
264 a sufficient number $N_L = (L/\delta\rho_f)^2$ of independent realizations of disorder, with
265 $\delta\rho_f \sim \lambda/(2NA)$ the standard diffraction-limited resolution. Indeed, the bias of
266 our aberration phase estimator is inversely proportional to $N_L^{26}$. In the present
267 case, the optimal value for $L$ has been found to be 18.6 $\mu$m.
268

269 The result of the aberration correction process is displayed in Fig. 4d at $z_t = 140$
270 $\mu$m. Strikingly, the estimated aberration laws exhibit strong phase fluctuations and
271 vary quickly between neighboring windows. This complex feature has two origins:
272 (i) the lateral fluctuations exhibited by the optical index inside the cornea; (ii)
273 the imperfections of the imaging system. The latter component accounts for the
274 difference observed between the input and output aberration transmittances (Sup-
275 plementary Section S4). In fact, the input aberration phase law accumulates not
276 only the input aberrations of the sample-arm but also those of the reference arm.
277 The sample-induced aberrations can be investigated independently from the im-
278 perfections of the experimental set up by considering the output aberration phase
279 matrix $\mathbf{\Phi}_{\text{out}}$ (Fig. 4d$_2$). Despite the prior compensation of the mean defocus by
280 minimizing the RPSF width, the aberration phase laws exhibit a residual defocus
281 that varies across the field-of-view due to the lateral variations of the optical in-
282 dex. Local shifts of the pupil function are also observed and result from a local
283 curvature of the coherence surface with respect to the focusing plane.

284 The extracted aberration phase laws can be used to estimate the transmission
285 matrices, $\mathbf{G}_{\text{in/out}}(z_t) = [G_{\text{in/out}}(\boldsymbol{\rho}_{\text{in/out}}, \boldsymbol{\rho}', z_t)]$, containing the impulse responses
286 $G_{\text{in/out}}(\boldsymbol{\rho}_{\text{in/out}}, \boldsymbol{\rho}', z_t)$ between the image voxels ($\boldsymbol{\rho}_{\text{in/out}}$) and the sample mapped
287 by the vector $\boldsymbol{\rho}'$ at each depth $z_t$ (Methods, Eq. 36). Aberrations are then com-
288 pensated by applying the phase conjugate of the transmission matrices³¹ at the
289 input and output of the focused **R**-matrix (Fig. 2a₃), such that:

$$
\mathbf{R}_{\rho\rho}''(z_t) = \mathbf{G}_{\text{out}}^{\dagger}(z_t) \times \mathbf{R}_{\rho\rho}(z_t) \times \mathbf{G}_{\text{in}}^{*}(z_t) \quad (8)
$$

290 The final image of the sample can be obtained by considering the diagonal elements
291 of the corrected matrix **R**''**ρρ**:

$$
I_M(\mathbf{r}_t) = R''(\boldsymbol{\rho}, \boldsymbol{\rho}, z_t). \qquad (9)
$$

292 Figures 2b3, c3 and d3 display longitudinal and transverse cross-sections of the
293 cornea obtained by RMI (see also Supplementary Movies 1, 2 and 3). The com-
294 parison with the confocal image [Figs. 2b2, c2 and d2] shows a clear gain in contrast
295 and resolution. This drastic improvement of the image can be quantified by exam-
296 ining the RPSF. While, at the previous step, the confocal peak exhibits a spreading
297 well beyond the diffraction limit and a background at depth due to forward multi-
298 ple scattering events (Figs. 2e1,f1), RMI compensates for these two issues and leads
299 to an almost ideal RPSF (Figs. 2e2,f2) over the whole field-of-view (Fig. 4e2). The
300 image improvement can be assessed by the RPSF (Fig. 2f3) whose Strehl ratio S
301 is increased by a factor 7 compared to its original value (Fig. 2f2). At the same
302 depth zt = 230 µm, the contrast χ is increased by a factor 11 and the resolution δρ
303 is enhanced by a factor four (Fig. 2f3), thus nearly reaching the diffraction-limited
304 value ρf ~ 500 nm.

305 The obtained three-dimensional image highlights several crucial features of the
306 cornea: its lamellar structure induced by the collagen fibrils (Fig. 2b3); (ii) the
307 complex network of nerves that covers the cornea; (iii) characteristic structures
308 of the cornea such as keratocytes and; (iv) stromal striae whose presence is an
309 indicator of keratoconus³². Such a high-resolution image can thus be of particular
310 importance for bio-medical diagnosis, given the high frame rate of our device.
Of course, RMI is not limited to the cornea but can be also applied to the deep
inspection of retina, skin or arteries, tissues whose structures are already monitored
by OCT but, until now, limited by a modest penetration depth.

Addressing Forward Multiple Scattering.

In that respect, Fig. 5 shows that RMI can also succeed in a stronger (forward)
315 multiple scattering regime, by imaging a resolution target through a more opaque
316 region of the cornea, closer to the iris (Fig. 5a, see Methods). Figures 5b and c
317 display the images obtained for a single illumination ($\mathbf{u}_{\text{in}} = \mathbf{0}$, Eq. 3) and after
318 the coherent compound of $N_{\text{in}} = 325$ illuminations (Eq. 5). While the holoscopic
319 image exhibits a random speckle due to multiple scattering (Fig. 5b), most of
320 the patterns of the resolution target are revealed by the digital confocal image
321 (Fig. 5c). However, it is poorly contrasted ($\chi \sim 0.9$ dB). This is confirmed by
322 the blurred feature of the RPSFs characteristic of multiple scattering, which high-
323 lights the cornea turbidity (Fig. 5g). The associated Strehl ratio is extremely low:
324 $S \sim 10^{-3}$. To tackle multiple scattering, a local analysis of wave distortions is
325 then performed (Methods, Eq. 33). The corresponding aberration phase laws are
326 displayed in Fig. 5e,f. They exhibit a complex feature characteristic of forward
327 multiple scattering with a broad spatial frequency content and a short-scale mem-
328 ory effect¹⁸. The phase conjugate of the associated transmission matrix (Eq. 8)
329 provides the final image (Eq. 9) displayed in Fig. 5d. The high contrast of the im-
330 age demonstrates the benefit of RMI. The comparison between original and final
331 RPSFs confirms the drastic improvement of the focusing quality (Figs. 5g,h). The
332 transverse resolution almost reaches the confocal resolution $\delta\rho_c$: $\delta\rho \sim 0.35$ µm.
333

The ability of RMI in overcoming high-order aberrations and multiple scatter-
ing thus constitutes a paradigm shift for deep optical microscopy. In its present
334 form, the penetration depth of RMI remains of the order of $\ell_t$, the typical depth
335 beyond which the memory effect vanishes³³. However, a multi-scale¹⁸ and/or multi-
336 conjugate³⁴ compensation of wave distortions can actually address even more com-
337
<figure><img src="image_6.png" /><figcaption>FIG. 5. **Imaging a resolution target through the cornea.** a, Experimental configuration. b, Holoscopy [scale bar: 20µm]. c, Digital confocal image. d, Matrix image based on the IPR algorithm. e, Maps of the aberration phase laws estimated from the pupil plane. f, Zoom on the white rectangle displayed in panel e. g, Maps of RPSFs obtained after the digital focusing process. h, Maps of RSPFs after aberration matrix compensation. Each RPSF is displayed over a de-scan area of 6.3 × 6.3 µm².</figcaption></figure>

plex scattering trajectories associated with extremely small isoplanatic patches (a
few µm). Moreover, the multi-spectral reflection matrix gives access to temporal
degrees of freedom that can be exploited for tailoring complex spatio-temporal
focusing laws²¹ required to overcome the diffusive limit. The mapping of the re-
fractive index will also be an important step to build accurate focusing laws inside
the medium³⁵. As shown by quantitative phase imaging of thin biological samples,
this physical parameter is also a quantitative marker for biology. Mapping the
refractive index in 3D and in an epi-detection geometry will pave the way towards
a quantitative imaging of biological tissues.

348

349 Discussion

350 In contrast with previous works that considered the reflection matrix at a single
351 frequency<sup>36</sup> or time-of-flight<sup>12,13,16</sup>, the measurement of a polychromatic reflection
352 matrix<sup>22</sup> allowed us to realize in post-processing: (i) a 3D confocal image of the
353 sample reflectivity on millimetric volumes (0.1 mm³ = 5 × 10⁹ voxels) in an ultra-
354 fast acquisition time (1.4 s) ; (ii) a local compensation of aberrations which usually
355 prevent deep imaging.

356 The required number of input wave-fronts depends on the aberration level and
357 scales as the number of resolution cells covered by the RPSF<sup>26</sup>. If this condi-
358 tion is not met, the diffuse halo extends beyond the RPSF extension (Fig. 5g)
359 and suppression of multiple scattering is then less efficient than it would be in
360 conventional confocal imaging (see the residual background in Fig. 5h). In prac-
361 tice, a compromise shall therefore be made between the frame rate and the image
362 quality. Note also that a great advantage of the matrix approach with respect to
363 conventional confocal imaging lies in its ability of local aberration compensation
364 in post-processing. Moreover, a spectral measurement and spatial multiplexing
365 of the wave-field provides a drastic enhancement of the signal-to-noise ratio by a
366 factor N<sub>ω</sub> × N<sub>in</sub>, i.e the number of spatio-temporal degrees of freedom provided
367 by our illumination sequence<sup>37</sup>.

368 An issue we have not considered in this paper is medium motion during the
369 acquisition of the reflection matrix. Of course, the assumption of a static medium is
370 everything but true especially for in-vivo applications<sup>38</sup>. To cope with the dynamic
371 features of the medium, two strategies can be followed. The first one is to limit the
372 measurement time of the **R**-matrix at its minimum, as allowed by our device using a
373 few illuminations. The second one is to develop algorithms that consider medium
374 motion during the measurement of **R**<sup>34</sup>. Interestingly, temporal fluctuations of
375 the medium's reflectivity and refractive index can provide a key information for
376 probing the multi-cellular dynamics in optical microscopy<sup>39,40</sup>.

377 Just as the concept of plane-wave imaging<sup>41</sup> revolutionized the field of ultra-
378 sound<sup>42,43</sup> by providing an unprecedented frame rate, our device will constitute
379 in a near future an ideal tool for probing the 3D dynamics of tissues at a much
380 smaller scale<sup>40,44</sup>. In medical imaging, an increase of the frame rate by one to two
381 orders of magnitude is necessary for ophtalmic applications due to eye movements
382 (10 Hz⁴⁵). It is also important in dermatology since standard OCT devices re-
383 main extremely slow to acquire 3D skin images (0.04 Hz⁴⁶). At last, our imaging
384 apparatus is an ideal tool for the monitoring of cell development in 3D whether
385 it be for embryology or pharmacology with the fast development of organoids⁴⁷.
386 For all these applications, the 3D imaging capability at a high frame rate will be
387 particularly decisive.
### 388 Methods

#### 389 Experimental components.

The following components were used in the experimental setup (Fig. 1): A swept laser source (800-875 nm; Superlum-850 HP), one galvanometer (Thorlabs, LSKGG4), one scan lens $L_1$ ($f_1$ = 110 mm), two immersion objective lenses (40×; NA, 0.8; Nikon), an imaging lens $L_2$ ($f_2$ = 250 mm) and an ultrafast camera (25 kHz; Phantom-v2640).

#### 395 Sample preparation.

In the presented experiment, the corneal sample was fixed with paraformaldehyde (4% concentration). Note that this substance tends to enhance scattering while absorption remains unchanged$^{48}$.

#### 399 Sampling of input and output wave-fields.

The dimension of the input pupil is $D_{in} \times D_{in} = 9 \times 9$ mm; the spatial sampling of input wave-fields is $\delta u_{in} = 600$ µm. Given the magnification of the output lens system (MO, $L_2$) system MO1 and the inter-pixel distance of the camera ($\delta s_{out} = 12$ µm), the output wave-field is sampled at a resolution close to $\lambda/(4NA)$, the theoretical achievable resolution for a confocal image$^{28}$: $\delta\rho_c = 290$ nm.

#### 405 Data acquisition and GPU processing.

All the interferograms of the acquisition sequence are recorded by the camera in 1.4 s and stored in its internal memory. Then, the whole data set (75 Gb) is transferred to the computer in 5 min. The numerical post-processing of the reflection matrix is performed by GPU (NVIDIA TITAN RTX) and takes 3.6 s per input wave-fronts. For the data set considered in this paper, all the focusing and aberration correction algorithms are performed in 1 hour.

#### 412 On-axis holography.
For each input wave-field, the interferogram recorded by the camera can be expressed as follows:

$$I(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = |E(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) + E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega)|^2 \quad (10)$$

with $E$ and $E_{\text{ref}}$, the wave-fields reflected by the sample and reference arms. Then a Fourier transform in the frequency domain is performed. The resulting intensity can be written as follows:

$$I(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, t) = E(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, t) \stackrel{t}{\circledast} E^*(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, -t) \quad (11)$$

$$+ E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, -t) \stackrel{t}{\circledast} E^*_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, -t) \quad (12)$$

$$+ E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, t) \stackrel{t}{\circledast} E^*(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, -t) \quad (13)$$

$$+ E(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, t) \stackrel{t}{\circledast} E^*_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, -t) \quad (14)$$

The two first terms (Eqs. 11 and 12) correspond to the self-interference of each arm with itself. Both contributions emerge at an optical depth close to zero ($t = 0$).

The two last terms correspond to the anti-causal (Eq. 13) and causal (Eq. 14) components of the interference between the two arms. By applying a Heavyside filter to $I(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, t)$ along the time dimension, one can isolate the causal contribution (Eq. 14). An inverse Fourier transform then yields the distorted wave-field:

$$D(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = E(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) E^*_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega). \quad (15)$$

If aberrations in the reference arm are neglected (Supplementary Section S3), the reference wave-field is a replica of the incident wave-field,

$$E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \exp \left( i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \mathbf{s}_{\text{out}} \right).$$
421 The multi-spectral reflection matrix is thus extracted using the following relation:

$$ \mathbf{R}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \mathbf{D}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \mathbf{s}_{\text{out}} \right). \quad (16) $$

422 Fresnel operators

423 The numerical focusing process is performed by means of Fresnel propagators.

424 For this purpose, the multi-spectral reflection matrix should be first projected in 425 the output pupil plane ($\mathbf{u}_{\text{out}}$) by a simple 2D spatial Fourier transform:

$$ \mathbf{R}_{uu}(\omega) = \mathbf{P}_{us}^*(\omega) \times \mathbf{R}_{su}(\omega) \quad (17) $$

426 where $\mathbf{P}_{us}(\omega) = [P(\mathbf{u}, \mathbf{s}, \omega)]$ is the Fourier transform operator:

$$ P(\mathbf{u}, \mathbf{s}, \omega) = e^{-j \frac{\omega \mathbf{u} \cdot \mathbf{s}}{c_0 f}} \quad (18) $$

427 with $f$, the focal length of the MOs and $c_0$ the vacuum light velocity. A Fresnel 428 phase law is then applied at the output of $\mathbf{R}_{uu}(\omega)$ to numerically shift the focal 429 plane, originally located in the middle of the sample ($z = 0$) to any depth $z$:

$$ \mathbf{R}_{\rho u}(z, \omega) = [\mathbf{P}_{u\rho}^{\top}(\omega) \circ \mathcal{F}_{u}(z, \omega)] \times \mathbf{R}_{uu}(\omega) \quad (19) $$

430 where the symbol $\circ$ accounts for the Hadamard (term-by-term) product. Each 431 column vector $\mathcal{F}_{u}(z, \omega)$ is a phase mask that accounts for the propagation of each 432 plane wave of transverse wave vector $\mathbf{k}_{\|} = \omega \mathbf{u}/(c_0 f)$ over a thickness $z$ of an 433 homogeneous medium of refractive index $n_0$:

$$ \mathcal{F}(\mathbf{u}, z, \omega) = e^{-j \left( \frac{n_0 \omega}{c_0} - k_z \right) z} \mathcal{O}(\mathbf{u}) \quad (20) $$
434 with

$$k_z = \frac{\omega}{c_0} \sqrt{n_0^2 - \frac{\|\mathbf{u}\|^2}{f^2}}, \qquad (21)$$

435 the longitudinal component of the wave vector, and $O(\mathbf{u})$, the finite pupil support: 436 $O(\mathbf{u}) = 1$ for $\|\mathbf{u}\| < fNA$ and zero elsewhere. Each reflection matrix $R_{\rho\mathbf{u}}(z, \omega) = 437 [R(\rho_{\text{out}}, \mathbf{u}_{\text{in}}, z, \omega)]$ connects each output virtual focusing point $\mathbf{r}_{\text{out}} = (\rho_{\text{out}}, z)$ to 438 each input illumination $\mathbf{u}_{\text{in}}$ at frequency $\omega$. The combination of Eqs. 17, 19 and 20 439 leads us to define the Fresnel operator $F_{\rho s}(z, \omega) = [F(\rho_{\text{out}}, \mathbf{s}_{\text{out}}, z, \omega)]$ that enables 440 the projection of the optical data from the camera sensors ($\mathbf{s}_{\text{out}}$) to any focal plane 441 $(\rho_{\text{out}})$ at expected depth $z$ (Eq. 1):

$$F_{\rho s}(z, \omega) = [\mathbf{P}_{\mathbf{u}\rho}^{\top}(\omega) \circ \mathcal{F}_{\mathbf{u}}(z, \omega)] \times \mathbf{P}_{\mathbf{us}}^{*}(\omega) \qquad (22)$$

442

The projection between the plane wave illumination basis ($\mathbf{u}_{\text{in}}$) and the focused 443 basis $(\rho_{\text{in}})$ can also be performed by means of a Fresnel propagator $F_{\rho u}$ this time 444 defined between the pupil plane and the each focal plane identified by their depth 445

446 $z$:

$$F_{\rho u}(z, \omega) = \mathbf{P}_{\mathbf{u}\rho}^{\dagger}(\omega) \circ \mathcal{F}_{\mathbf{u}}(z, \omega) \qquad (23)$$

447

448 Local estimation of focusing quality

To probe the local RPSF, the field-of-view is divided at each effective depth 449 $z_t$ into regions that are defined by their central midpoint $\rho_p = (x_p, y_p)$ and their 450 spatial extension $L$. A local average of the back-scattered intensity can then be 451 performed in each region:

$$RPSF_l(\Delta\rho, \rho_p, z, z_t) = \langle |R(\rho + \Delta\rho/2, \rho - \Delta\rho/2, z, t)|^2 W_L(\rho - \rho_p) \rangle_{\rho} \qquad (24)$$
453 where the symbol $\langle \cdots \rangle_m$ stands for an average over the variable $m$ in subscript.
454 $W_L(\rho - \rho_p) = 1$ for $|x - x_p| < L/2$ and $|y - y_p| < L/2$, and zero otherwise. In
455 this paper, a spatial window of size $L = 18.6$ µm has been used to smooth out
456 fluctuations due to the sample inhomogeneous reflectivity$^{29}$.

457 Strehl ratio and contrast

458 On the one hand, a quantification of the focusing quality can be peformed with
459 the RPSF by computing the associated Strehl ratio $S$. In reflection, this quantity
460 can actually be defined as the ratio between the confocal energy and the overall
461 backscattered energy, such that:

$$
S = \frac{RSPSF(\Delta\rho = 0)}{\sum_{\Delta\rho} RSPSF(\Delta\rho)}. \qquad (25)
$$

462 On the other hand, the image contrast $\chi$ can be evaluated by considering the ratio
463 between the confocal peak of the RPSF and its diffuse background:

$$
\chi = \frac{RPSF(\Delta\rho = 0)}{\langle RSPSF(\Delta\rho) \rangle_{||\Delta\rho|| > 2\delta\rho}}. \quad (26)
$$

464 where the diffuse background is here estimated by averaging the RPSF outside the
465 confocal peak, *i.e* arbitrarily beyond $2\delta\rho$.

466 Local compensation of wave-distortions

467 The starting point is the time-gated reflection matrix $\mathbf{R}_{\rho\rho}(z_t)$, obtained after
468 tuning the focusing plane and coherence volume at each echo time $t$. The first
469 step is a projection of $\mathbf{R}_{\rho\rho}(z_t)$ in the pupil plane at input via a numerical Fourier
470 transform:

$$
\mathbf{R}_{\rho u}(z_t) = \mathbf{R}_{\rho \rho}(z_t) \times \mathbf{P}_{u \rho}^{\top}(\omega_c) \quad (27)
$$

471 with $\omega_c$, the central frequency. An input distortion matrix is then built by per-
472 forming a element-wise product between $\mathbf{R}_{u\rho}(z_t)$ and the phase conjugate reference
$$ \mathbf{D}_{\rho u}(z_t) = \mathbf{R}_{\rho u}(z_t) \circ \mathbf{P}_{\mathbf{u}\rho}^{\dagger}(\omega_c) \quad (28) $$

$$ C_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{u}'_{\text{in}}, \mathbf{r}_{\mathbf{p}}) = \langle D(\rho_{\text{out}}, \mathbf{u}_{\text{in}}, z_t) D(\rho_{\text{out}}, \mathbf{u}'_{\text{in}}, z_t) W_L(\rho_{\text{out}} - \rho_{\mathbf{p}}) \rangle_{\rho_{\text{out}}} \quad (29) $$

$$ \mathbf{R}'_{\rho\rho}(z_t) = \{\mathbf{R}_{\rho u}(z_t) \circ \exp[-i\Phi_{\text{in}}(z_t)]\} \times \mathbf{P}_{\mathbf{u}\rho}^* \quad (30) $$

$$ \mathbf{R}'_{\mathbf{u}\rho}(z_t) = \mathbf{P}_{\mathbf{u}\rho} \times \mathbf{R}'_{\rho\rho}(z_t). \quad (31) $$

$$ \mathbf{D}_{\mathbf{u}\rho}(z_t) = \mathbf{P}_{\mathbf{u}\rho}^* \circ \mathbf{R}'_{\mathbf{u}\rho}(z_t) \quad (32) $$

$$ C_{\text{out}}(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \mathbf{r}_{\mathbf{p}}) = \langle D(\mathbf{u}_{\text{out}}, \rho_{\text{in}}, z_t) D_{\text{out}}^*(\mathbf{u}'_{\text{out}}, \rho_{\text{in}}, z_t) W_L(\rho_{\text{in}} - \rho_{\mathbf{p}}) \rangle_{\rho_{\text{in}}} \quad (33) $$

<footnote n="473">473 matrix $\mathbf{P}_{\mathbf{u}\rho}(\omega_c)$ that would be obtained in absence of aberrations¹⁵ (Supplementary Section S6):</footnote>

<footnote n="475">475 A local correlation matrix $\mathbf{C}_{\text{in}}$ of wave distortions is then built around each point 476 $\mathbf{r}_{\mathbf{p}} = (\rho_{\mathbf{p}}, z_t)$ of the field-of-view (Supplementary Section S7). Its coefficients write:</footnote>

<footnote n="477">477 Iterative phase reversal (see further) is then applied to each correlation matrix 478 $\mathbf{C}_{\text{in}}(\mathbf{r}_{\mathbf{p}})^{18}$ (Supplementary Section S8). The resulting input phase laws, $\Phi_{\text{in}}(z_t) = 479 [\phi_{\text{in}}(\mathbf{u}_{\text{in}}, \rho_{\mathbf{p}}, z_t)]$, are used to compensate for the wave distortions undergone by the 480 incident wave-fronts:</footnote>

<footnote n="481">481 The corrected matrix $\mathbf{R}'_{\rho\rho}$ is only intermediate since phase distortions undergone 482 by the reflected wave-fronts remain to be corrected.</footnote>

<footnote n="483">483 To that aim, $\mathbf{R}'_{\rho\rho}(z_t)$ is now projected in the pupil plane at output:</footnote>

<footnote n="484">484 An output distortion matrix is then built:</footnote>

<footnote n="485">485 From $\mathbf{D}_{\mathbf{u}\rho}$, one can build a correlation matrix $\mathbf{C}_{\text{out}}$ for each point $\mathbf{r}_{\mathbf{p}}$:</footnote>
486 The IPR algorithm described further is then applied to each matrix $C_{\text{out}}(\mathbf{r}_{\text{p}})$. 487 The resulting output phase laws, $\mathbf{\Phi}_{\text{out}}(z_t) = [\phi_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{p}}, z_t)]$, are leveraged to 488 compensate for the residual wave distortions undergone by the reflected wave- 489 fronts:

$$ \mathbf{R}_{\rho\rho}''(z_t) = \mathbf{P}_{u\rho}^{\dagger} \times \{\exp[-i\mathbf{\Phi}_{\text{out}}(z_t)] \circ \mathbf{R}_{u\rho}'(z_t)\} \quad (34) $$

490 Iterative phase reversal algorithm.

491 The IPR algorithm is a computational process that provides an estimator of 492 the phase of the transmittance that links each point $\mathbf{u}$ of the pupil plane with each 493 voxel $\mathbf{r}_{\text{p}} = (\boldsymbol{\rho}_{\text{p}}, z_t)$ of the cornea volume$^{18}$. To that aim, the correlation matrix 494 $\mathbf{C}$ computed over the spatial window $W_L$ centered around a given point $(\boldsymbol{\rho}_{\text{p}}, z_t)$ 495 is considered (Eqs. 29 and 33). Mathematically, the algorithm is based on the 496 following recursive relation:

$$ \mathbf{\Phi}_{\text{in/out}}^{(n)}(\boldsymbol{\rho}_{\text{p}}, z_t) = \arg \left\{ \mathbf{C}_{\text{in/out}}(\boldsymbol{\rho}_{\text{p}}, z_t) \times \exp \left[ i \mathbf{\Phi}_{\text{in/out}}^{(n-1)}(\boldsymbol{\rho}_{\text{p}}, z_t) \right] \right\} \quad (35) $$

497 where $\mathbf{\Phi}_{\text{in/out}}^{(n)}$ is the estimator of the transmittance phase at the $n^{\text{th}}$ iteration of 498 the phase reversal process. $\mathbf{\Phi}_{\text{in/out}}^{(0)}$ is an arbitrary wave-front that initiates the 499 process (typically a flat phase law) and $\phi_{\text{in/out}} = \lim_{n\to\infty} \phi_{\text{in/out}}^{(n)}$ is the result of 500 IPR.

501

502 Transmission matrices.

503 The transmission matrices $\mathbf{G}_{\text{in/out}}$ used to deconvolve the focused $\mathbf{R}$-matrix 504 (Eq. 8) can be deduced from the estimated aberration phase laws $\mathbf{\Phi}_{\text{in/out}}$ as follows 505 (Eqs. 30 and 34):

$$ \mathbf{G}_{\text{in/out}}(z_t) = \{\mathbf{P}_{u\rho}^T \circ \exp[i\mathbf{\Phi}_{\text{in/out}}(z_t)]\} \times \mathbf{P}_{u\rho}^*. \quad (36) $$

506
508 **Resolution target experiment.**

509

510 A resolution target is placed behind the cornea in the focal plane of the micro-
511 scope objective (NA=0.8). The reflection matrix **R** is measured with N_in = 325
512 plane waves, corresponding to a full scan of the immersion MO pupil (NA=0.8,
513 refractive index n₀ = 1.33). The interferograms are recorded by N_out = 512² pixels
514 of the camera, corresponding to an output FOV of Ω_out × Ω_out = 179.2 × 179.2
515 µm², with a spatial sampling δρ_out = 350 nm. Finally, N_ω = 180 independent
516 frequencies are used to probe the sample. A RPSF maximization is performed to
517 tune the coherence plane with the focal plane in order to get the digital confocal
518 image displayed in Fig. 5c. A local compensation of wave distortions is performed
519 over reduced spatial windows of size L = 9.3 µm in order to get the final image
520 displayed in Fig. 5d.

521

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

524

525 **Code availability.** Codes used to post-process the optical data within this paper
526 are available upon request at Zenodo (https://zenodo.org/records/10674114)

527

528 **Acknowledgments.** The authors wish to thank A. Badon for initial discussions
529 on the project, K. Irsch for providing the corneal sample, A. Le Ber for providing
530 the iterative phase reversal algorithm and F. Bureau for his help on Supplementary
531 movies.

532

533 **Funding Information.** The authors are grateful for the funding provided by
534 the European Research Council (ERC) under the European Union's Horizon 2020
535 research and innovation program (grant agreement nos. 610110 and 819261, 536 HELMHOLTZ* and REMINISCENCE projects, respectively). This project has 537 also received funding from Labex WIFI (Laboratory of Excellence within the 538 French Program Investments for the Future; ANR-10-LABX-24 and ANR-10- 539 IDEX-0001-02 PSL*) and from CNRS Innovation (Prematuration program, MA- 540 TRISCOPE project).

541

542 **Author Contributions.** A.A. initiated and supervised the project. P.B., V.B. 543 and A.A. designed the experimental setup. P.B., V.B. and N.R. built the experi- 544 mental set up. V.B. designed the acquisition scheme. V.B., P.B., N.G. and E.A. 545 developed the post-processing tools. P.B. performed the corneal imaging experi- 546 ment. V.B., P.B. and A.A. analyzed the experimental results. V.B., P.B. and A.A. 547 performed the theoretical study. P.B., V.B. and A.A. prepared the manuscript. 548 P.B., V.B., N.G., A.C.B., M.F. and A.A. discussed the results and contributed to 549 finalizing the manuscript.

550

551 **Competing interests.** P.B., V.B., M.F., C.B. and A.A. are named inventors on 552 french patent FR2207334 (filing date 18.07.2022), which is related to the techniques 553 described in this Article.
# Supplementary Material

This document provides further information on: (i) the experimental set up; (ii) 3D images of the cornea; (iii) the theoretical expression of the multi-spectral reflection matrix; (iv) the theoretical expression of the focused reflection matrix; (v) the reflection point spread function; (vi) the local distortion matrix; (vii) the corresponding correlation matrix; (viii) iterative phase reversal; (ix) the proof-of-concept experiment with a resolution target placed behind the cornea.

## 561 S1. EXPERIMENTAL SET UP

The full experimental set up is displayed in Fig. S1. Compared with Fig. 1 of the accompanying paper, it shows the presence of a scan lens and of a tube lens in order to focus the incident light in the pupil plane of the microscope objectives at input. It also highlights the control of light polarization in order to minimize spurious reflections. The beam splitter is polarized and quarter wave plates are placed in the two arms such that the reflected light in the two arms is transmitted to the CMOS camera. The amount of light injected in the two arms is controlled by means of two polarizers, P1 and P2, placed before the polarized beam splitter. An analyzer A1 allows us to project the sample and reference beams on the same polarization and make them interfere in the camera plane.

## 572 S2. OTHER 3D IMAGES OF THE CORNEA

In complement of Fig. 2 of the accompanying paper, Figs. S2b₁-d₁ shows longitudinal and transverse cross-sections of the cornea obtained via spectra domain OCT (Fig. S2a₁). Its comparison with the holoscopic image displayed in Figs. 2b₁-d₁ of the accompanying paper illustrate the effect of numerical focusing. While the
<figure><img src="image_7.png" /><figcaption>FIG. S1. Detailed experimental set up. P: Polarizers; A: Analyzer; QW: quarter wave plates.</figcaption></figure>

577 OCT image is completely blurred due to multiple scattering and finite depth-of-field of the high-NA microscope objective, the numerical back-focusing of the 578 optical wave-field improves the image contrast. Nevertheless, although the bright-579 est are revealed, the holoscopic image still suffers from a strong multiple scatter-580 ing background that generates a random speckle. This speckle can be smoothed 581 by an incoherent average of the holoscopic image obtained for each illumination 582 (Figs. S2a₂-d₂). However, such an incoherent compound image remains largeley 583
584 perfectible since the contrast remains very weak. On the contrary, a coherent 585 compound of holoscopic images provides a much contrasted image of the cornea 586 (Figs. 2a₂-d₂ of the accompanying paper). Indeed, a coherent combination of 587 multi-illuminations acts as a confocal pinhole.

### 588 S3. MULTI-SPECTRAL REFLECTION MATRIX

589 In this section, we express theoretically the multi-spectral reflection matrix 590 recorded by the experimental set up of Fig. S1. To that aim, we will rely on a 591 simple Fourier optics model to describe the multi-spectral reflection matrix. For 592 the sake of simplicity, this model is scalar.

593 The wave field $E_s(s_{\text{out}}, \mathbf{u}_{\text{in}}, \omega)$ reflected by the sample arm in the camera plane 594 can be expressed as follows:

$$E_s(s_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = S(\omega) \iiint G(s_{\text{out}}, \boldsymbol{\rho}_s, z_s, \omega) \gamma(\boldsymbol{\rho}_s, z_s) E_{\text{in}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s, \omega) d\boldsymbol{\rho}_s dz_s. \quad (S1)$$

595 $S(\omega)$ is the amplitude of light source at frequency $\omega$. $G(s_{\text{out}}, \boldsymbol{\rho}_s, z_s)$, the Green's 596 function between the sample mapped by the vector $(\boldsymbol{\rho}_s, z_s)$ and the CCD sensors 597 identified by $s_{\text{out}}$. $\gamma(\boldsymbol{\rho}_s, z_s)$ represents the sample reflectivity. $E_{\text{in}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z, \omega)$ 598 describes the incident wave-field that can be expressed as follows:

$$E_{\text{in}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s, \omega) = S(\omega) \mathcal{F}(\mathbf{u}_{\text{in}}, z_s, \omega) T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z) \exp \left[ -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_s \right] \quad (S2)$$

599 where $\mathcal{F}(\mathbf{u}_{\text{in}}, z_s, \omega)$ is the Fresnel phase law that describes plane wave propagation 600 inside an homogeneous medium of effective index $n$ (Eq. 20 in the accompanying paper). The transmission matrix $T = [T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z)]$ accounts for the wave 601 distortions induced by the fluctuations of the optical index inside the medium.

602 In the reference arm, a mirror is placed in the focal plane of the MO and displays 603 a uniform reflectivity: $\gamma(\boldsymbol{\rho}_m, z_m) = \gamma_m \delta(z_m)$, with $\delta$ the Dirac distribution and $\gamma_m$ 604
<figure><img src="image_8.png" /><figcaption>FIG. S2. Intermediate 3D images of the cornea. **a**, Imaging Methods. **b**, B-scan image showing one longitudinal section of the cornea reflectivity. **c-d**, En-face image of the cornea at *z* = 150 µm and 275 µm, respectively [scale bar: 75 µm]. In panels (a)-(d), subscripts 1 and 2 stand for spectral domain OCT and multi-illumination holoscopy (incoherent compound).</figcaption></figure>
the mirror surface reflectivity. The reference wave-field is thus given by:

$$E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \gamma_m S(\omega) \iint G_{\text{ref}}(\mathbf{s}_{\text{out}}, \boldsymbol{\rho}_m) T_{\text{ref}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_m) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_m \right) d\boldsymbol{\rho}_m. \quad (\text{S3})$$

with $G_{\text{ref}}(\mathbf{s}_{\text{out}}, \boldsymbol{\rho}_m, \omega)$, the Green's functions between the focal plane of the MO ($\boldsymbol{\rho}_m$) and the CCD sensors ($\mathbf{s}_{\text{out}}$) and $T_{\text{ref}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_m)$, the transfer function describing the aberrations undergone by the incident wave in the reference arm due to experimental imperfections (MO, misalignment, etc.). Assuming isoplanicity in the reference arm, the Green's function $G_{\text{ref}}$ can be replaced by a spatially-invariant impulse response $H_{\text{ref}}$ between the focal plane and the CCD sensors: $G_{\text{ref}}(\mathbf{s}_{\text{out}}, \boldsymbol{\rho}_m) = H_{\text{ref}}(\boldsymbol{\rho}_m + \mathbf{s}_{\text{out}})$. Under the same hypothesis, the transfer function $T_{\text{ref}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_m)$ becomes an aberration transmittance $\mathcal{T}_{\text{ref}}(\mathbf{u}_{\text{in}})$, defined as the Fourier transform of the reference arm point spread function $H_{\text{ref}}$:

$$\mathcal{T}_{\text{ref}}(\mathbf{u}_{\text{in}}) = \iint d\boldsymbol{\rho} H_{\text{ref}}(\boldsymbol{\rho}) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho} \right). \quad (\text{S4})$$

Under the isoplanatic assumption, Equation S3 thus simplifies into:

$$E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \gamma_m S(\omega) \mathcal{T}_{\text{ref}}(-\mathbf{u}_{\text{in}}) \mathcal{T}_{\text{ref}}(\mathbf{u}_{\text{in}}) \exp \left( i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \mathbf{s}_{\text{out}} \right) \quad (\text{S5})$$

If aberrations in the reference arm are neglected, we retrieve the fact the reference wave-field is a replica of the input wave-front:

$$E_{\text{ref}}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \gamma_m S(\omega) \exp \left( i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \mathbf{s}_{\text{out}} \right). \quad (\text{S6})$$

In the following, we will not make this assumption and will consider the more general expression of $E_{\text{ref}}$ given in Eq. S5.

The coefficients of the multi-spectral matrix $\mathbf{R}_{\text{su}}(\omega)$ are recorded by isolating the interference between the sample beam, $E_s$, and the reference beam, $E_{\text{ref}}$
622 (Eqs. 15 and 16 of the accompanying paper):

$$R(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = E(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) E_{\text{ref}}^{*}(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) \exp \left( i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \mathbf{s}_{\text{out}} \right) \quad (S7)$$

623 Using Eqs. S1, S2, S5, the last equation can be rewritten as follows:

$$R(\mathbf{s}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \gamma_m |S(\omega)|^2 \iiint G(\mathbf{s}_{\text{out}}, \boldsymbol{\rho}_s, z_s, \omega) \gamma(\boldsymbol{\rho}_s, z_s) \\ \times T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s) \mathcal{F}(\mathbf{u}_{\text{in}}, z_s, \omega) \\ \times \mathcal{T}_{\text{ref}}^{*}(-\mathbf{u}_{\text{in}}) \mathcal{T}_{\text{ref}}^{*}(\mathbf{u}_{\text{in}}) \exp \left[ -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_s \right] d\boldsymbol{\rho}_s dz_s \quad (S8)$$

624 One can already notice from this expression that the aberrations undergone by 625 the reference wave-field ($\mathcal{T}_{\text{ref}}^{*}(-\mathbf{u}_{\text{in}})\mathcal{T}_{\text{ref}}^{*}(\mathbf{u}_{\text{in}})$) emerge at the input of the recorded 626 reflection matrix.

627 S4. FOCUSED REFLECTION MATRIX

628 In this section, we describe theoretically the numerical focusing process leading 629 to a time-gated focused reflection matrix at each depth of the sample.

630 First, a spatial Fourier transform over the camera pixels $\mathbf{s}_{\text{out}}$ leads to the reflec- 631 tion matrix $\mathbf{R}_{\text{uu}}$ in the pupil basis (Eqs. 17 and 18 of the accompanying paper):

$$R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) = \gamma_m |S(\omega)|^2 \iiint T(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z) \mathcal{F}(\mathbf{u}_{\text{out}}, z_s, \omega) \gamma(\boldsymbol{\rho}_s, z_s) \\ \times T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z) \mathcal{F}(\mathbf{u}_{\text{in}}, z_s, \omega) \\ \times \mathcal{T}_{\text{ref}}^{*}(-\mathbf{u}_{\text{in}}) \mathcal{T}_{\text{ref}}^{*}(\mathbf{u}_{\text{in}}) \quad (S9) \\ \times \exp \left[ -i \frac{2\pi}{\lambda f} (\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}}) \cdot \boldsymbol{\rho}_s \right] d\boldsymbol{\rho}_s dz_s. \quad (S10)$$

632 As for incident light (Eq. S2), the return path is decomposed in the plane wave 633 basis as the product between a Fresnel phase law $\mathcal{F}(\mathbf{u}_{\text{out}}, z_s, \omega)$, accounting for
634 free-space wave propagation in an homogeneous medium of refractive index $n$, and
635 the transfer function $T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s)$ that grasps the wave distortions induced by the
636 refractive index fluctuations such that:

$$ T(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z_s) \mathcal{F}(\mathbf{u}_{\text{out}}, z_s, \omega) = \iint G(\mathbf{s}_{\text{out}}, \boldsymbol{\rho}_s, z_s) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{out}} \cdot \mathbf{s}_{\text{out}} \right) d\boldsymbol{\rho}_s. \quad (S11) $$

637 Numerical focusing at depth $z$ (Eqs. 1 and 4 of the accompanying paper) then
638 consists in compensating wave diffraction by applying the phase conjugate of the
639 Fresnel propagator for refractive index $n_0$ at input and output before a spectral
640 Fourier transform:

$$ R(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z, t) = \gamma_m \sum_{\mathbf{u}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} \int d\omega |S(\omega)|^2 \exp(i\omega t) \mathcal{F}_0^*(\mathbf{u}_{\text{out}}, z, \omega) R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{in}}, z, \omega) \\ \times \exp \left[ i \frac{2\pi}{\lambda f} (\mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_{\text{in}} + \mathbf{u}_{\text{out}} \cdot \boldsymbol{\rho}_{\text{out}}) \cdot \boldsymbol{\rho}_s \right] $$

Injecting Eq. S9 leads to the following expression for the coefficients of $\mathbf{R}_{\rho\rho}(z, \omega)$:

$$ \begin{aligned} R(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z, t) = & \gamma_m \sum_{\mathbf{u}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} \int d\omega |S(\omega)|^2 \exp(i\omega t) \\ & \times \iiint T(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z) \gamma(\boldsymbol{\rho}_s, z_s, \omega) T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s) \mathcal{T}_{\text{ref}}^*(-\mathbf{u}_{\text{in}}) \mathcal{T}_{\text{ref}}^*(\mathbf{u}_{\text{in}}) \\ & \times \mathcal{F}(\mathbf{u}_{\text{out}}, z_s, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{out}}, z_s, \omega) \mathcal{F}(\mathbf{u}_{\text{in}}, z_s, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{in}}, z_s, \omega) \\ & \times \exp \left\{ -i \frac{2\pi}{\lambda f} \left[ \mathbf{u}_{\text{in}} \cdot (\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in}}) + \mathbf{u}_{\text{out}} \cdot (\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{out}}) \right] \right\} d\boldsymbol{\rho}_s dz_s. \end{aligned} \quad (S12) $$

641 The positions of the coherence volume and focusing plane are determined by the
642 cancellation of the Fresnel phase laws,

$$
\begin{equation}
\begin{split}
\mathcal{F}(\mathbf{u}_{\text{in}}, z, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{in}}, z_s, \omega) \mathcal{F}(\mathbf{u}_{\text{out}}, z, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{out}}, z_s, \omega) ={}& \\
& \exp \left[ i \frac{\omega}{c_0} \left( \sqrt{n_0^2 - \frac{\|\mathbf{u}_{\text{in}}\|^2}{f^2}} + \sqrt{n_0^2 - \frac{\|\mathbf{u}_{\text{out}}\|^2}{f^2}} - 2n_0 \right) z \right] \\
& \times \exp \left[ -i \frac{\omega}{c_0} \left( \sqrt{n^2 - \frac{\|\mathbf{u}_{\text{in}}\|^2}{f^2}} + \sqrt{n^2 - \frac{\|\mathbf{u}_{\text{out}}\|^2}{f^2}} \right) z_s \right].
\end{split}
\tag{S13}
\end{equation}
$$

643 Under the paraxial approximation, these Fresnel phase laws can be developed as
644 follows:

$$
\begin{equation}
\begin{split}
\mathcal{F}(\mathbf{u}_{\text{in}}, z, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{in}}, z_s, \omega) \mathcal{F}(\mathbf{u}_{\text{out}}, z, \omega) \mathcal{F}_0^*(\mathbf{u}_{\text{out}}, z_s, \omega) ={}& \\
& \exp \left[ -2i \frac{n\omega}{c_0} z_s \right] \exp \left[ -i \frac{\omega}{c_0} \left( \frac{\|\mathbf{u}_{\text{in}}\|^2}{2f^2} + \frac{\|\mathbf{u}_{\text{out}}\|^2}{2f^2} \right) \left( \frac{z}{n_0} - \frac{z_s}{n} \right) \right],
\end{split}
\tag{S14}
\end{equation}
$$

The cancellation of the first phase term defines the real position of the coherence
volume $z_s = z_t = c_0t/(2n)$ that appears at an effective depth $z_0 = c_0t/(2n_0) =$
$(n/n_0)z_t$ (Fig. S3a). Previous expression of $R(\rho_{\text{out}}, \rho_{\text{in}}, z, t)$ (Eq. S12) can be
rewritten as follows:

$$
R(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z, t) = \gamma_m f(t - 2nz_t/c) \stackrel{t}{\otimes} \quad (\text{S15})
$$

$$
\begin{equation}
\begin{split}
& \sum_{\omega} \sum_{\mathbf{u}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} \iiint T(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z) \gamma(\boldsymbol{\rho}_s, z_s) \\
& \qquad \times T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_s, z_s) \mathcal{T}_{\text{ref}}^{*}(-\mathbf{u}_{\text{in}}) \mathcal{T}_{\text{ref}}^{*}(\mathbf{u}_{\text{in}}) \exp(i\omega t) \\
& \qquad \times \exp \left[ -i \frac{\omega}{c_0} \left( \frac{\|\mathbf{u}_{\text{in}}\|^2}{2f^2} + \frac{\|\mathbf{u}_{\text{out}}\|^2}{2f^2} \right) \left( \frac{z}{n_0} - \frac{z_s}{n} \right) \right] \\
& \qquad \times \exp \left\{ -i \frac{2\pi}{\lambda f} \left[ \mathbf{u}_{\text{in}} \cdot (\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in}}) + \mathbf{u}_{\text{out}} \cdot (\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{out}}) \right] \right\} d\boldsymbol{\rho}_s dz_s.
\end{split}
\tag{S16}
\end{equation}
$$

645 with $f(t) = \int d\omega |S(\omega)|^2 e^{i\omega t}$, the time response of the microscope. The symbol $\otimes$
646 stands for convolution over variable t. The position $z_f$ of the focusing plane is 647 obtained when the parabolic phase term cancels in the previous expression, that 648 is to say for $z_s = z_f = (n_0/n)z_t$. The apparent defocus induced by the mismatch 649 between n and $n_0$ is thus equal to:

$$ \Delta z = z_f - z_0 = z_t(n_0/n - n/n_0) \qquad (S18) $$

650 An index mismatch thus implies a defocus distance $\Delta z$ that increases linearly with 651 $z_t$ (Fig. S3a). In the present study, the estimated defocus is roughly constant with 652 $z_t$. It thus means that the cornea displays an effective optical index $n \sim n_0 = 1.33$ 653 and that the observed defocus rather originates from the different lengths between 654 the sample and reference arms.

Once the focusing plane is matched with the coherence volume, the Fresnel phase laws in Eq. S15 vanish. Assuming $n = n_0$, the coefficients of the time-gated focused reflection matrix $R_{\rho\rho}(z_t)$ can be derived as follows:

$$ R(\rho_{\text{out}}, \rho_{\text{in}}, z_t) = \gamma_m \sum_{\mathbf{u}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} \iint T_{\text{out}}(\mathbf{u}_{\text{out}}, \rho_s, z_t) \gamma(\rho_s, z_t) T_{\text{in}}(\mathbf{u}_{\text{in}}, \rho_s, z_t) \\ \times \exp \left\{ i \frac{2\pi}{\lambda f} [\mathbf{u}_{\text{in}} \cdot (\rho_s - \rho_{\text{in}}) + \mathbf{u}_{\text{out}} \cdot (\rho_s - \rho_{\text{out}})] \right\} d\rho_s. \quad (S19) $$

with

$$ T_{\text{out}}(\mathbf{u}_{\text{out}}, \rho_s, z_t) \equiv T(\mathbf{u}_{\text{out}}, \rho_s, z_t) $$

and

$$ T_{\text{in}}(\mathbf{u}_{\text{in}}, \rho_s, z_t) \equiv T(\mathbf{u}_{\text{in}}, \rho_s, z_t) T_{\text{ref}}^{*}(-\mathbf{u}_{\text{in}}) T_{\text{ref}}^{*}(\mathbf{u}_{\text{in}}). $$

655 While the output transmission matrix $T_{\text{out}}$ corresponds to the sample transmission 656 matrix $T$, the input transmission matrix $T_{\text{in}}$ grasps both the sample and reference 657 arm aberrations.

658 The Fourier transform of the transmission coefficients $T_{\text{in/out}}$ in Eq. S19 provide
<figure><img src="image_9.png" /><figcaption>FIG. S3. Mismatch between the coherence volume and focusing plane. a. For a medium of refractive index $n > n_0$, the focusing plane at $z_f$ is shifted from the coherence volume at $z_t$ and expected ballistic depth $z_0$. b. A defocus $\Delta z$ can be applied in post-processing in order to make coincide the coherence and focusing planes.</figcaption></figure>

659 local PSFs, $H_{in/out}$, such that:

$$H_{\text{in/out}}(\boldsymbol{\rho}, \boldsymbol{\rho}_s, z_t) = \sum_{\mathbf{u}} T_{\text{in/out}}(\mathbf{u}, \boldsymbol{\rho}_s, z_t) \exp\left(i \frac{2\pi}{\lambda f} \mathbf{u} \cdot \boldsymbol{\rho}\right). \quad (\text{S20})$$

660 The time-gated reflection matrix can be rewritten as follows:

$$R(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z_t) = \iint d\boldsymbol{\rho}_s H_{\text{out}}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_s, z_t) \gamma(\boldsymbol{\rho}_s, z) H_{\text{in}}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_s, z_t) \quad (\text{S21})$$

661 The latter expression can be recast as a function the impulse responses $G_{\text{in/out}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_{\text{in/out}}, z_t)$

662 between input/output focusing points and points $\boldsymbol{\rho}_s$ mapping the sample. Both

663 quantities are actually linked as follows:

$$G_{\text{in/out}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_{\text{in/out}}, z_t) = H_{\text{in/out}}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in/out}}, \boldsymbol{\rho}_s, z_t) \quad (\text{S22})$$
664 Injecting the last expression into Eq. S21 leads to the following expression:

$$R(\rho_{\text{out}}, \rho_{\text{in}}, z_t) = \iint G_{\text{out}}(\rho_{\text{out}}, \rho_s, z_t) \gamma(\rho_s, z_t) G_{\text{in}}(\rho_{\text{in}}, \rho_s, z_t) \quad (\text{S23})$$

665 Under a matrix formalism, the last expression can be rewritten as follows:

$$\mathbf{R}_{\rho\rho}(z_t) = \mathbf{G}_{\text{out}}(z_t) \times \mathbf{\Gamma}(z_t) \times \mathbf{G}_{\text{in}}^{\top}(z_t), \quad (\text{S24})$$

666 $\Gamma$ describes the scattering process inside the sample. Under a single scattering 667 assumption, this matrix is diagonal. Its coefficients then correspond the sample 668 reflectivity $\gamma(\rho_s, z_t)$. $\mathbf{G}_{\text{in}}$ and $\mathbf{G}_{\text{out}}$ are the input and output focusing matrices. 669 Their coefficients, $G_{\text{in/out}}(\rho_{\text{in/out}}, \rho_s, z_t)$, describe the transverse amplitude distribution of the focal spot when trying to focus at point $(\rho_{\text{in/out}}, z_t)$. 670

## 671 S5. REFLECTION POINT SPREAD FUNCTION

672 As mentioned in the accompanying paper, the off-diagonal coefficients of $\mathbf{R}_{\rho\rho}(z, t)$ 673 enable to probe the focusing quality at any voxel by investigating the reflection 674 point spread function (RPSF, Eq. 6). To express theoretically the latter quantity, a 675 local isoplanatic assumption shall be made. This hypothesis implies that the PSFs 676 $H_{\text{in/out}}$ are locally invariant by translation. This leads us to define local spatially- 677 invariant PSFs $h_{\text{in/out}}$ around each central midpoint $\rho_p$ at each time-of-flight such 678 that:

$$H_{\text{in/out}}(\rho_s - \rho_{\text{in/out}}, \rho_s, z, t) = h_{\text{in/out}}(\rho_s - \rho_{\text{in/out}}, \rho_p, z_p, t). \quad (\text{S25})$$

679 The second assumption is to consider the medium reflectivity $\gamma(\rho_s, z)$ as random:

$$\langle \gamma(\rho_1, z) \gamma^*(\rho_2, z) \rangle = \langle |\gamma|^2 \rangle \delta(\rho_2 - \rho_1), \quad (\text{S26})$$
680 By combining those assumptions with Eq. S21, the ensemble average of RPSF(Δρ, ρp, z, t)

681 (Eq. 7) can be expressed as follows:

$$ \langle \text{RPSF}(\Delta\boldsymbol{\rho}, \boldsymbol{\rho}_p, z, t) \rangle = \langle |\gamma|^2 \rangle \times \left[ |h_{\text{in}}|^2 \stackrel{\Delta\boldsymbol{\rho}}{\otimes} |h_{\text{out}}|^2 \right] (\Delta\boldsymbol{\rho}, \boldsymbol{\rho}_p, z_p, t). \quad (\text{S27}) $$

## 682 S6. LOCAL DISTORTION MATRIX

683 Wave distortions can be investigated both at input and output of the reflection matrix (see Methods of the accompanying paper). Here we will consider the properties of the output distortion matrix but the same theoretical developments can be made at input.

687 The output distortion matrix can be built by first projecting the time-gated focused reflection matrix $R_{\rho\rho}(z_t)$ in the pupil plane at output:

$$ \mathbf{R}_{u\rho}(z_t) = \mathbf{P}_{u\rho} \times \mathbf{R}_{\rho\rho}(z_t). \quad (\text{S28}) $$

689 Then, the distorted component of the wave-field can be extracted by subtracting the geometric phase expected in an ideal case (without aberrations). Mathematically, this can be performed using the following matrix element wise product:

$$ \mathbf{D}_{\text{out}}(z_t) = \mathbf{R}_{u\rho}(z_t) \circ \mathbf{P}_{u\rho}^*. \quad (\text{S29}) $$

692 or in terms of matrix coefficients,

$$ D(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z_t) = \sum_{\boldsymbol{\rho}_{\text{out}}} R(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z_t) \exp \left[ -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{out}} \cdot (\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_{\text{in}}) \right]. \quad (\text{S30}) $$
Injecting Eqs. S20 and S21 into the last equation yields

$$D(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z_t) = \iint T_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z_t) \gamma(\boldsymbol{\rho}_s, z_t) H_{\text{in}}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_s, z) \quad (S31)$$

$$\times \exp \left[ -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{out}} \cdot (\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in}}) \right] d\boldsymbol{\rho}_s. \quad (S32)$$

In previous papers$^{15,18}$, we showed that the distortion matrix **D** highlights spatial correlations of the reflected wave-field induced by the shift-shift memory effect$^{30,33}$: Waves produced by nearby points inside an anisotropic scattering medium generate highly correlated distorted wave-fields in the pupil plane. 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.

To extract and exploit this local memory effect for imaging, the field-of-illumination should be subdivided into overlapping regions$^{18}$ that are defined by their central midpoint ($\boldsymbol{\rho}_p, z_t$) and their lateral extension $L$. All of the distorted components associated with focusing points $\boldsymbol{\rho}_{\text{in}}$ located within each region are extracted and stored in a local distortion matrix $D'_{\text{out}}(\boldsymbol{\rho}_p, z_t)$:

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

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

Under a local isoplanatic assumption (Eq. S25), the aberrations can be modelled by a local transmittance $T_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_p, z_t)$ around each point ($\boldsymbol{\rho}_p, z_t$), such that $T_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_s, z_t) \simeq T_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_p, z_t)$. This transmittance is the Fourier transform of the local PSF $h_{\text{out}}(\boldsymbol{\rho}, \boldsymbol{\rho}_p, z_t)$:

$$T_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_p, z_t) = \iint h_{\text{out}}(\boldsymbol{\rho}, \boldsymbol{\rho}_p, z) \exp \left( \frac{2\pi}{\lambda f} \mathbf{u}_{\text{out}} \cdot \boldsymbol{\rho} \right) d\boldsymbol{\rho} \quad (S34)$$
Under this assumption, Eq.S31 can be rewritten as follows:

$$D'(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{p}}, z_t) = \underbrace{\mathcal{T}_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{p}}, z_t)}_{\text{transmittance}} \times \underbrace{\iint \gamma(\boldsymbol{\rho} + \boldsymbol{\rho}_{\text{out}}, z_t)h_{\text{in}}(\boldsymbol{\rho}, \boldsymbol{\rho}_{\text{p}}, z_t) \exp\left(-i\frac{2\pi}{\lambda f}\mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}\right) d\boldsymbol{\rho}}_{\text{virtual source}} \quad (S35)$$

709 The physical meaning of this last equation is the following: Each distorted wave- 710 field corresponds to the diffraction of a virtual source synthesized inside the 711 medium modulated by the transmittance $\mathcal{T}_{\text{out}}$ of the sample between the focal 712 and pupil planes. Each virtual source is spatially incoherent due to the random 713 reflectivity of the medium, and its size is governed by the spatial extension of 714 the input focal spot. The idea is now to smartly combine each virtual source to 715 generate a coherent guide star and estimate $\mathcal{T}_{\text{out}}$ independently from the sample 716 reflectivity.

## 717 S7. CORRELATION OF WAVE DISTORTIONS

718 To do so, the correlation matrix $\mathbf{C}_{\text{out}} = \mathbf{D}_{\text{out}}\mathbf{D}_{\text{out}}^{\dagger}$ is an excellent tool. Its 719 coefficients write as follows

$$C_{\text{out}}(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_{\text{p}}, z_t) = N_{\text{W}}^{-1} \sum_{\boldsymbol{\rho}_{\text{in}}} D'(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{p}}, z_t) D'^{*}(\mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{p}}, z_t) \quad (S36)$$

720 The matrix $\mathbf{C}_{\text{out}}(\boldsymbol{\rho}_{\text{p}}, z_t)$ can be decomposed as the sum of its ensemble average, 721 the covariance matrix $\langle \mathbf{C}_{\text{out}} \rangle (\boldsymbol{\rho}_{\text{p}}, z_t)$, and a perturbation term $\delta\mathbf{C}_{\text{out}}(\boldsymbol{\rho}_{\text{p}}, z_t)$:

$$\mathbf{C}_{\text{out}}(\boldsymbol{\rho}_{\text{p}}, z_t) = \langle \mathbf{C}_{\text{out}} \rangle (\boldsymbol{\rho}_{\text{p}}, z_t) + \delta\mathbf{C}_{\text{out}}(\boldsymbol{\rho}_{\text{p}}, z_t). \quad (S37)$$
722 The intensity of the perturbation term scales as the inverse of the number $N_L =$ 723 $(L/\delta\rho_0)^2$ of resolution cells in each sub-region$^{49,50}$:

$$ \langle |\delta C_{\text{out}}(\mathbf{u}, \mathbf{u}', \boldsymbol{\rho}_p, z_t)|^2 \rangle = \frac{\langle |C_{\text{out}}(\mathbf{u}, \mathbf{u}, \boldsymbol{\rho}_p, z_t)|^2 \rangle}{N_L} \quad (\text{S38}) $$

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

726 Under assumptions of local isoplanicity (Eqs. S25 and S35) and random reflec- 727 tivity,

$$ \langle \gamma(\boldsymbol{\rho}_s, z) \gamma^*(\boldsymbol{\rho}'_s, z) \rangle = \langle |\gamma|^2 \rangle \delta(\boldsymbol{\rho}_s - \boldsymbol{\rho}'_s). \quad (\text{S39}) $$

728 The coefficients of the covariance matrix can be expressed as follows$^{51}$:

$$ \langle \mathbf{C}_{\text{out}} \rangle (\boldsymbol{\rho}_p, z_t) = [\mathcal{T}_{\text{out}}(\boldsymbol{\rho}_p, z_t) \circ \mathbf{P}_{\mathbf{u}\boldsymbol{\rho}}] \times \mathbf{C}_H(\mathbf{r}_p) \times [\mathbf{P}_{\mathbf{u}\boldsymbol{\rho}} \circ \mathcal{T}_{\text{out}}(\boldsymbol{\rho}_p, z_t)]^{\dagger}, \quad (\text{S40}) $$

or in terms of matrix coefficients,

$$ \begin{aligned} \langle \mathbf{C}_{\text{out}} \rangle (\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_p, z_t) &= \mathcal{T}_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_p, z_t) \mathcal{T}^*_{\text{out}}(\mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_p, z_t) \\ & \quad \times \iint d\boldsymbol{\rho}_s |h_{\text{in}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_p, z_t)|^2 \exp \left[ -i \frac{2\pi}{\lambda f} (\mathbf{u}_{\text{out}} - \mathbf{u}'_{\text{out}}) \cdot \boldsymbol{\rho}_s \right] \\ &= \mathcal{T}_{\text{out}}(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_p, z_t) \mathcal{T}^*_{\text{out}}(\mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_p, z_t) \\ & \quad \times \underbrace{\left[ \mathcal{T}_{\text{in}} \mathbf{u}^* \mathcal{T}_{\text{in}} \right] (\mathbf{u}_{\text{out}} - \mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_p, z_t)}_{=C_H(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \boldsymbol{\rho}_p, z_t)} \end{aligned} \quad (\text{S41}) $$

729 where the symbol * stands for correlation product over variable $\mathbf{u}$. $C_H$ is a refer- 730 ence correlation matrix that would be measured in an homogeneous cornea for a 731 virtual reflector whose scattering distribution corresponds to the output focal spot 732 intensity $|h_{\text{in}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_p, z_t)|^2$. The covariance matrix $\langle \mathbf{C}_{\text{in}} \rangle (\boldsymbol{\rho}_p, z_t)$ thus corresponds 733 to the same experimental situation but for a virtual reflector embedded into the 734 heterogeneous cornea under study.
735 S8. ITERATIVE PHASE REVERSAL

736 An estimator $\exp[i\Phi_{\text{out}}(\rho_p, z_t)]$ of the local aberration transmittance $T_{\text{out}}(\rho_p, z_t)$ 737 can be extracted by applying an iterative phase reversal algorithm to $C_{\text{out}}$ (Eq. 34 738 of the accompanying paper). It mimics an iterative time reversal process on the 739 virtual reflector described above but imposes a constant amplitude for the time- 740 reversal invariant. This iterative phase reversal (IPR) process converges towards a 741 wavefront that maximizes the coherence of the wave-field reflected by the virtual 742 reflector¹⁸.

743 This IPR process assumes the convergence of the correlation matrix $C_{\text{out}}$ 744 (Eq. S37) towards its ensemble average $\langle C_{\text{out}} \rangle$, the covariance matrix⁵⁰,⁵¹. In fact, 745 this convergence is never fully realized and $C_{\text{out}}$ should be decomposed as the sum 746 of this covariance matrix $\langle C_{\text{out}} \rangle$ and the perturbation term $\delta C_{\text{out}}$ (Eq. S37). This 747 incomplete convergence towards the covariance matrix leads to a phase error $\delta\phi$ 748 made on our estimation of each aberration phase law. The variance of this error 749 scales as follows:

$$ \langle |\delta\phi_{\text{out}}(\mathbf{u})|^2 \rangle \sim \frac{1}{N_L C_{\text{in}}} \qquad (\text{S42}) $$

750 with $C_{\text{in}}$, the coherence factor that is a direct indicator of the focusing quality⁵², 751 ranging from 0 for a fully blurred guide star to 2/9 for a diffraction-limited focal 752 spot⁵³. On the one hand, this scaling of the phase error with $N_W$ explains why 753 each spatial window should be large enough to encompass a sufficient number of 754 independent realizations of disorder. On the other hand, it should be small enough 755 to grasp the spatial variations of aberration phasr laws across the field-of-view. A 756 compromise has thus to be found between these two opposite requirements. It has 757 led us to to choose spatial windows of size $L = 18.6 \ \mu\text{m}$ to ensure the convergence of 758 the IPR process. Note that, contrary to a recent study performed in an extremely 759 opaque cornea¹⁸, a multi-scale approach of wave distortions is not necessary here 760 since the cornea under study displays smoother fluctuations of optical index over
larger characteristic length scales. As previous works dealing with the CLASS algorithm<sup>14,16</sup>, these imaging conditions ensure:

• A lower level of transverse aberrations, hence a higher coherence factor that allows a direct convergence of IPR over reduced spatial windows.

• Larger isoplanatic patches, hence a direct convergence of the IPR process for reduced spatial windows and no need to iterate the IPR process over smaller spatial windows.

## 768 S9. IMAGING A RESOLUTION TARGET THROUGH AN OPAQUE REGION OF THE CORNEA

In addition to the three-dimensional imaging experiment of the cornea presented in the accompanying manuscript, an academic experiment has been performed and consists in imaging a resolution target behind the same cornea. This experiment provides a validation of the whole RMI process with a ground-truth object. It also allows us to test RMI in a stronger scattering regime since we image the resolution target through a more opaque region of the cornea, closer to the iris. This way, we demonstrate the ability of RMI in overcoming high-order aberrations and forward multiple scattering.

Figure S4a shows the experimental configuration in which a resolution target is placed behind the cornea in the focal plane of the microscope objective (Methods). Figure S3b displays the image obtained for a single illumination ($u_{in} = 0$, Eq. 3). This image exhibits a random speckle due to the multiple scattering events induced by the heterogeneities of the cornea. A direct combination of each backscattered wave-field for each plane wave illumination yields a raw confocal image (Fig. S4c, Eq. 5). Despite spatial multiplexing and time-gating, this image does not provide a better result because of the mismatch between the coherence volume and the focal plane (Fig. S3a). An optimization of the RPSF allows us to finely tune these
<figure><img src="image_10.png" /><figcaption>FIG. S4. Imaging a resolution target through the cornea. a, Experimental configuration. b, Holoscopy [scale bar: 20µm]. c, Digital confocal image. d, Matrix image based on a SVD approach¹⁵. e, Matrix image based on the IPR algorithm.</figcaption></figure>

two planes (Fig. S3b), thereby providing the digital confocal image displayed in Fig. S4d. This image reveals most of the patterns of the resolution target. However, it is poorly contrasted ($\chi \sim 0.9$ dB) because of the aberrations induced by the fluctuations of the refractive index inside the cornea. The corresponding RPSFs display a blurred feature characteristic of multiple scattering, which highlights the cornea turbidity (Fig. 5g). The associated Strehl ratio is extremely low: $S \sim 10^{-3}$

Figure S4e shows a corrected image based on a SVD decomposition of the full-field **D**-matrix as proposed in a previous paper¹⁵. This image is built from the 25th first eigenstates of **D**. It displays a better contrast than the optimized confocal image (Fig. S4d) but it remains far from being perfect due to the high number of isoplanatic patches in the FOV. Note that considering more eigenstates would not improve the image contrast because the higher-rank eigenstates are polluted
799 by multiple scattering noise.

800 To circumvent this problem, an analysis of wave distortions can be performed
801 over reduced spatial windows (Methods). The corresponding aberration phase
802 laws are displayed in Fig. 5e,f. They exhibit a complex feature characteristic of
803 forward multiple scattering with a broad spatial frequency content and a short-
804 scale memory effect¹⁸. The phase conjugate of the associated transmission matrix
805 (Eq. 8) provides the final image (Eq. 9) displayed in Fig. S4f. The high contrast of
806 the image demonstrates the benefit of a local analysis of wave distortions compared
807 with a global SVD approach (Fig. S4e). The comparison between original and
808 final RPSFs confirms the drastic improvement of the focusing quality (Figs. 5g
809 and h). The Strehl ratio is increased by a factor 14. Nevertheless, the final
810 RPSFs are not perfect: (i) The contrast remains limited ($\chi \sim 10$); (ii) The
811 transverse resolution is slightly larger than the confocal resolution $\delta\rho_c$: $\delta\rho \sim 0.35$
812 $\mu$m. Aberrations associated with an isoplanatic length smaller than $L$ cannot
813 be addressed by RMI and could explain, at least partially, the residual flaws of
814 the final RPSFs. Moreover, the incomplete sampling of the illumination sequence
815 truncates the diffuse halo as displayed by the original RPSFs (Fig. 5g). This lack of
816 information alters the aberration correction process and accounts for the residual
817 background displayed by the final RPSFs (Fig. 5b). This academic experiment
818 thus highlights both the benefits and limits of RMI.

<footnote n="819">819 [1] V. Ntziachristos, Going deeper than microscopy: The optical imaging frontier in biology, Nat. Methods 7, 603 (2010).</footnote>

<footnote n="821">821 [2] S. Gigan, Optical microscopy aims deep, Nat. Photonics 11, 14 (2017).</footnote>

<footnote n="822">822 [3] J. Bertolotti and O. Katz, Imaging in complex media, Nat. Phys. 18, 1008 (2022).</footnote>

<footnote n="823">823 [4] M. J. Booth, Adaptive optical microscopy: the ongoing quest for a perfect image, Light Sci. Appl. 3, e165 (2014).</footnote>
[5] T.-W. Wu and M. Cui, Numerical study of multi-conjugate large area wavefront correction for deep tissue microscopy, Opt. Express **23**, 7463 (2015).

[6] J.-H. Park, Z. Yu, K. Lee, P. Lai, and Y. Park, Perspective: Wavefront shaping techniques for controlling multiple light scattering in biological tissues, APL Photonics **3**, 100901 (2018).

[7] T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, Interferometric synthetic aperture microscopy, Nat. Phys. **3**, 129 (2007).

[8] S. G. Adie, B. W. Graf, A. Ahmad, P. S. Carney, and S. A. Boppart, Computational adaptive optics for broadband optical interferometric tomography of biological tissue, Proc. Nat. Acad. Sci. USA **109**, 7175 (2012).

[9] A. Ahmad, N. D. Shemonski, S. G. Adie, H.-S. Kim, W.-M. W. Hwu, P. S. Carney, and S. A. Boppart, Real-time in vivo computed optical interferometric tomography, Nat. Photonics **7**, 444 (2013).

[10] D. Hillmann, H. Spahr, C. Hain, H. Sudkamp, G. Franke, C. Pfäffle, C. Winter, and G. Hüttmann, Aberration-free volumetric high-speed imaging of in vivo retina, Sci. Rep. **6**, 35209 (2016).

[11] O. Haim, J. Boger-Lombard, and O. Katz, Image-guided computational holographic wavefront shaping, arXiv: 2305.12232 (2023).

[12] S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J.-S. Lee, Y.-S. Lim, Q.-H. Park, and W. Choi, Imaging deep within a scattering medium using collective accumulation of single-scattered waves, Nat. Photonics **9**, 253 (2015).

[13] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, Smart optical coherence tomography for ultra-deep imaging through highly scattering media, Sci. Adv. **2**, e1600370 (2016).

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

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

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

[17] Y. Kwon, J. H. Hong, S. Kang, H. Lee, Y. Jo, K. H. Kim, S. Yoon, and W. Choi, Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin, Nat. Commun. **14**, 105 (2023).

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

[19] D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, C. Puliafito, et al., Optical coherence tomography, Science **254**, 1178 (1991).

[20] A. Badon, A. C. Boccara, G. Lerosey, M. Fink, and A. Aubry, Multiple scattering limit in optical microscopy, Opt. Express **25**, 28914 (2017).

[21] Y.-R. Lee, D.-Y. Kim, Y. Jo, M. Kim, and W. Choi, Exploiting volumetric wave correlation for enhanced depth imaging in scattering medium, Nat. Commun. **14** (2023).

[22] Y. Zhang, D. Minh, Z. Wang, T. Zhang, T. Chen, and C. W. Hsu, Deep imaging inside scattering media through virtual spatiotemporal wavefrontshaping, arXiv:2306.08793 (2023).

[23] B. Povazay, A. Unterhuber, B. Hermann, H. Sattmann, H. Arthaber, and W. Drexler, Full-field time-encoded frequency-domain optical coherence tomography, Opt. Express **14**, 7661 (2006).
[24] D. Hillmann, G. Franke, C. Lührs, P. Koch, and G. Hüttmann, Efficient holoscopy image reconstruction, Opt. Express **20**, 21247 (2012).

[25] E. Auksorius, D. Borycki, P. Stremplewski, K. Lizewski, S. Tomczewski, P. Niedzwiedziuk, B. L. Sikorski, and M. Wojtkowski, In vivo imaging of the human cornea with high-speed and high-resolution fourier-domain full-field optical coherence tomography, Biomed. Opt. Exp. **11**, 2849 (2020).

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

[27] P. Stremplewski, E. Auksorius, P. Wnuk, L. Kozon, P. Garstecki, and M. Wojtkowski, In vivo volumetric imaging by crosstalk-free full-field OCT, Optica **6**, 608 (2019).

[28] A. Sentenac and J. Mertz, Unified description of three-dimensional optical diffraction microscopy: from transmission microscopy to optical coherence tomography: tutorial, J. Opt. Soc. Am. A **35**, 748 (2018).

[29] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry, Reflection Matrix Approach for Quantitative Imaging of Scattering Media, Phys. Rev. X **10**, 021048 (2020).

[30] B. Judkewitz, R. Horstmeyer, I. M. Vellekoop, I. N. Papadopoulos, and C. Yang, Translation correlations in anisotropically scattering media, Nat. Phys. **11**, 684 (2015).

[31] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media, Phys. Rev. Lett. **104**, 100601 (2010).

[32] K. Grieve, D. Ghoubay, C. Georgeon, G. Latour, A. Nahas, K. Plamann, C. Crotti, R. Bocheux, M. Borderie, T.-M. Nguyen, F. Andreiuolo, M.-C. Schanne-Klein, and V. Borderie, Stromal striae: a new insight into corneal physiology and mechanics, Sci. Rep. **7**, 13584 (2017).
[33] G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, Generalized optical memory effect, Optica 4, 886 (2017).

[34] M. Kang, W. Choi, W. Choi, and Y. Choi, Fourier holographic endoscopy for imaging continuously moving objects, Opt. Express 31, 11705 (2023).

[35] M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, Multi-layer born multiple-scattering model for 3d phase microscopy, Optica 7, 394 (2020).

[36] S. M. Popoff, A. Aubry, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, Exploiting the Time-Reversal Operator for Adaptive Optics, Selective Focusing, and Scattering Pattern Analysis, Phys. Rev. Lett. 107, 263901 (2011).

[37] J. Mertz, Introduction to optical microscopy (Cambridge University Press, Cambridge, United Kingdom; New York, NY, 2019) Chap. 12, p. 254, second edition ed.

[38] M. Jang, H. Ruan, I. M. Vellekoop, B. Judkewitz, E. Chung, and C. Yang, Relation between speckle decorrelation and optical phase conjugation (OPC)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skin, Biomed. Opt. Express 6, 72 (2014).

[39] C. Apelian, F. Harms, O. Thouvenin, and A. C. Boccara, Dynamic full field optical coherence tomography: subcellular metabolic contrast revealed in tissues by interferometric signals temporal analysis, Biomed. Opt. Express 7, 1511 (2016).

[40] J. Scholler, K. Groux, O. Goureau, J.-A. Sahel, M. Fink, S. Reichman, C. Boccara, and K. Grieve, Dynamic full-field optical coherence tomography: 3d live-imaging of retinal organoids, Light Sci. Appl. 9, 140 (2020).

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

[42] E. Macé, G. Montaldo, I. Cohen, M. Baulac, M. Fink, and M. Tanter, Functional ultrasound imaging of the brain, Nat. Methods 8, 662 (2011).
933 [43] C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, and M. Tan-
ter, Ultrafast ultrasound localization microscopy for deep super-resolution vascular
imaging, Nature **527**, 499 (2015).

936 [44] M. Munter, M. vom Endt, M. Pieper, M. Casper, M. Ahrens, T. Kohlfaerber, 937 R. Rahmanzadeh, P. Konig, G. Huttmann, and H. Schulz-Hildebrandt, Dynamic 938 contrast in scanning microscopic oct, Opt. Lett. **45**, 4766 (2020).

939 [45] V. Mazlin, P. Xiao, K. Irsch, J. Scholler, K. Groux, K. Grieve, M. Fink, and A. C. 940 Boccara, Optical phase modulation by natural eye movements: application to time- 941 domain FF-OCT image retrieval, Biomed. Opt. Exp. **13**, 902 (2022).

942 [46] F. Latriglia, J. Ogien, C. Tavernier, S. Fischman, M. Suppa, J.-L. Perrot, and 943 A. Dubois, Line-field confocal optical coherence tomography (LC-OCT) for skin 944 imaging in dermatology, Life **13**, 2268 (2023).

945 [47] S. P. Harrison, S. F. Baumgarten, R. Verma, O. Lunov, A. Dejneka, and G. J. 946 Sullivan, Liver organoids: Recent developments, limitations and potential, Front. 947 Med. **8**, 574047 (2021).

948 [48] M. F. G. Wood, N. Vurgun, M. A. Wallenburg, and I. A. Vitkin, Effects of formalin 949 fixation on tissue optical polarization properties, Phys. Med. Biol. **56**, N115 (2011).

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

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

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

959 [52] R. Mallart and M. Fink, Adaptive focusing in scattering media through sound?speed
960

inhomogeneities: The van Cittert Zernike approach and focusing criterion, J.

961

Acoust. Soc. Am. **96**, 3721 (1994).

962

[53] S. Silverstein, Ultrasound scattering model: 2-d cross-correlation and focusing

963

criteria-theory, simulations, and experiments, IEEE Trans. Ultrason. Ferroelectr.

964

Freq. Control **48**, 1023 (2001)

965