domain: phys spi

text:
Paul Balondrade

Victor Barolle

Nicolas Guigui

Emeric Auriant

Nathan Rougier

Claude Boccara

Mathias Fink

Alexandre Aubry
email: alexandre.aubry@espci.fr

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

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 ℓ s (∼100 µm in biological tissues).

To circumvent this issue, adaptive optics (AO) has been transposed from astronomy to microscopy for the last twenty years [START_REF] Booth | Adaptive optical microscopy: the ongoing quest for a perfect image[END_REF] . 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 ℓ t (∼ 1 mm in biological tissues). Multi-conjugate AO could increase the corrected field-ofview [START_REF] Wu | Numerical study of multi-conjugate large area wavefront correction for deep tissue microscopy[END_REF] , but this would be at the price of a much more complex optical setup and an extremely long optimization process [START_REF] Park | Perspective: Wavefront shaping techniques for controlling multiple light scattering in biological tissues[END_REF] .

More recently, following seminal works that proposed post-processing computational strategies for AO [START_REF] Ralston | Interferometric synthetic aperture microscopy[END_REF][START_REF] Adie | Computational adaptive optics for broadband optical interferometric tomography of biological tissue[END_REF][START_REF] Ahmad | Real-time in vivo computed optical interferometric tomography[END_REF][START_REF] Hillmann | Aberration-free volumetric high-speed imaging of in vivo retina[END_REF][START_REF] Haim | Image-guided computational holographic wavefront shaping[END_REF] , a reflection matrix approach has been developed for deep imaging [START_REF] Kang | Imaging deep within a scattering medium using collective accumulation of single-scattered waves[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF][START_REF] Kwon | Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin[END_REF][START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] . 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 diffractionlimited resolution for each pixel of the field-of-view. Nevertheless, the existing approaches suffer from several limitations. In most experimental works [START_REF] Kang | Imaging deep within a scattering medium using collective accumulation of single-scattered waves[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF][START_REF] Kwon | Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin[END_REF][START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] , the reflection matrix is time-gated around the ballistic time as usually performed in time-domain OCT [START_REF] Huang | Optical coherence tomography[END_REF] . Such a measurement has one main advantage since it enables the temporal filtering of most multiply-scattered photons [START_REF] Badon | Multiple scattering limit in optical microscopy[END_REF] . However, it also suffers from two strong drawbacks. First, time-gating means that a large part of the information on the medium is discarded: Only the weakly distorted paths are recorded and can be compensated by a spatial phase modulation of the incident and reflected wave-fronts. Second, volumetric imaging can only be obtained by a mechanical axial scanning of the sample, which limits the frame rate F ps to, at best, 10 6 voxels.s -1 for a high quality correction over millimetric FOVs.

To go beyond, an acquisition of a spectral reflection matrix is required in order to capture all the information required for the three-dimensional imaging of a sample. In recent works [START_REF] Lee | Exploiting volumetric wave correlation for enhanced depth imaging in scattering medium[END_REF][START_REF] Zhang | Deep imaging inside scattering media through virtual spatiotemporal wavefrontshaping[END_REF] , the spatio-temporal degrees of freedom exhibited by the reflection matrix have been exploited for tailoring dispersive focusing laws.

However, the acquisition rate was slow (F ps ∼ 10 3 voxels.s -1 ) because the number of input wave-fronts scaled as the number of voxels in the image. Moreover, the experimental demonstration was limited to the imaging of a resolution target through a scattering medium [START_REF] Lee | Exploiting volumetric wave correlation for enhanced depth imaging in scattering medium[END_REF][START_REF] Zhang | Deep imaging inside scattering media through virtual spatiotemporal wavefrontshaping[END_REF] or a sparse medium made of colloidal particles [START_REF] Zhang | Deep imaging inside scattering media through virtual spatiotemporal wavefrontshaping[END_REF] .

In this paper, we go beyond an academic proof-of-concept and address the extremely challenging case of ultra-fast 3D imaging of biological tissues themselves (nerves, cells, collagen, extracellular matrix etc.). In particular, we will show how the number of input wave-fronts can be drastically decreased by deterministic focusing operations applied to the reflection matrix guided by a self-portrait of the focusing process.

To that aim, we report on a measurement of the multi-spectral reflection matrix at a much higher frame rate (F ps ∼ 10 9 voxels.s -1 ), with a 3D imaging demonstration on an ex-vivo opaque cornea at a resolution of 0.29 µm and 0.5 µm in the transverse and axial directions, respectively. The experimental set up combines a Fourier-domain full-field OCT (FD-FF-OCT) setup [START_REF] Považay | Full-field time-encoded frequency-domain optical coherence tomography[END_REF][START_REF] Hillmann | Efficient holoscopy image reconstruction[END_REF][START_REF] Auksorius | In vivo imaging of the human cornea with high-speed and high-resolution fourier-domain full-field optical coherence tomography[END_REF] with a coherent 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 spectral Fourier transform and numerical refocusing can provide a 3D image of the sample for each incident wave-front [START_REF] Považay | Full-field time-encoded frequency-domain optical coherence tomography[END_REF][START_REF] Hillmann | Efficient holoscopy image reconstruction[END_REF][START_REF] Auksorius | In vivo imaging of the human cornea with high-speed and high-resolution fourier-domain full-field optical coherence tomography[END_REF] 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 aberrations. Interestingly, reflection matrix imaging (RMI) can go beyond by decoupling input and output focusing points at each time-of-flight. A focused reflection matrix is synthesized and measures the cross-talk between each point inside the sample.

While previous works only considered focusing points at the same depth [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF] , we show here that their axial scan gives access to a self-portrait of the light focusing process. A minimization of the point spread function extension enables an autofocus process at each depth of the sample. Finally, a compensation of transverse aberrations is performed by means of a local analysis of wave distortions [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] . A digital clearing of long-scale refractive index 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.

Results

Recording the Multi-Spectral Reflection Matrix.

3D matrix imaging is based on the measurement of a multi-spectral reflection matrix from the scattering sample. The experimental setup and procedure are described in Fig. 1 (see Methods and Supplementary Figure S1). Inspired by spectral domain FFOCT [START_REF] Stremplewski | In vivo volumetric imaging by crosstalk-free full-field OCT[END_REF] , it simply consists in a Linnik interferometer (Fig. 1a). In the first arm, a reference mirror is placed in the focal plane of a microscope objective (MO). The second arm contains the scattering sample to be imaged through an identical MO. This interferometer is illuminated by a swept source through two scanning mirrors and a lens that allows a raster scanning of the focal spot in the MO pupil planes (Fig. 1a 1 ). The sample and reference mirror are thus illuminated by a set of plane waves at each frequency of the light source bandwidth (Fig. 1a 2 ).

The reflected waves are collected through the same MOs and, ultimately, interfere on a camera conjugated with the focal plane. For each input wave-front of coordinate u in in the pupil plane, the interferogram I(s out , u in , ω) recorded at frequency ω (Fig. 1a 3 ) provides one column of the reflection matrix R su = [R(s out , u in , ω)] (see Methods and Fig. 1b), where s out is the transverse location of each camera sensor.

In the opaque cornea experiment, the reflection matrix R is measured with N in = 177 plane waves, corresponding to a full scan of the immersion MO pupil (NA=0.8, refractive index n 0 = 1.33). The interferograms are recorded by N out = 1024 2 pixels of the camera, corresponding to an output FOV of Ω out ×Ω out = 297× 297 µm 2 , with a spatial sampling δρ out = 290 nm. Finally, N ω = 201 independent frequencies are used to probe the sample within the frequency bandwidth [800; 875] nm of the light source. All the information about the sample is thus contained in the 10 10 coefficients acquired in 1.4 s. In the following, we show how to postprocess this wealth of optical data to build a 3D highly-contrasted image of the cornea at a diffraction-limited resolution.

Ultra-fast Three-Dimensional Imaging.

To that aim, the most direct path is to perform, a Fourier transform over frequency ω of the back-scattered wave-field recorded for one illumination 23 : This is the principle of FF-SS-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 [START_REF] Barolle | Manifestation of aberrations in full-field optical coherence tomography[END_REF] , which is prohibitory for 3D imaging. 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 [START_REF] Hillmann | Aberration-free volumetric high-speed imaging of in vivo retina[END_REF] . This numerical focusing process is performed by means of Fresnel propagators.

(δz f ∼ n 0 λ/{2 sin 2 [asin(N A/2)]} ∼ 2.7 µm)
For this purpose, the multi-spectral reflection matrix is projected at output in the focused basis:

R ρu (z, ω) = F * us (z, ω) × R su (ω) (1)
where the symbol * stands for phase conjugate. F us (z, ω) = [F (u, s, z, ω)] is the Fresnel operator that describes free-space propagation from the camera (s) to any focal plane (ρ) located at expected depth z in the sample (Methods). Each frequency component of R ρu (z, ω) should then be recombined in order to time gate the singly-scattered photons. In practice, an inverse Fourier transform over frequency ω is performed and yields an R-matrix as a function of photon's timeof-flight t:

R ρu (z, t) = dωR ρu (z, ω)e jωt . (2)
At each time t, the single scattering contribution of the wave-field corresponds to photons that have been scattered in a coherence volume located at a depth z t in the sample and of thickness δz t ∼ c 0 λ 2 /(2n 0 ∆λ) ∼ 3.5 µm. When the focusing plane and the coherence volume coincide (Fig. 2a 1 ), an holoscopic image of the sample, I H , can be obtained for each input wave-front u in (Fig. 2a 1 ):

I H (r t , u in ) = R(ρ out , u in , z t , t). (3)
with r t = (ρ out , z t ). In practice, an exact matching between the focusing plane and coherence volume is difficult to obtain especially for deep imaging (i.e low single-to-multiple scattering ratio). We will describe further how matrix imaging can provide a robust observable for this fine tuning. S2). To get rid of it, the single-to-multiple scattering ratio shall be increased [START_REF] Badon | Multiple scattering limit in optical microscopy[END_REF] . For this purpose, a spatial filtering of multiply-scattered photons can be performed by means of a confocal filter. Nevertheless, this operation is extremely sensitive to the focusing quality inside the sample. A prior optimization of the focusing process is thus needed.

Digital confocal imaging.

To that aim, the dual reflection matrix is projected in the focused basis both at input and output. Mathematically, it simply consists in a numerical input focusing of R ρu using the Fresnel propagator F uρ that describes free space propagation from the MO pupil plane and the focal plane at depth z (Methods):

R ρρ (z, t) = dωR ρu (z, ω) × F † uρ (z, ω)e -jωt . (4)
where the symbol † stands for transpose conjugate. Expressed in the focused basis, the reflection matrix R ρρ (z, t) contains the responses at each time-of-flight t between virtual sensors of expected positions r in = (ρ in , z) and r out = (ρ out , z).

The focused R-matrix is equivalent to the time-gated reflection matrix considered in previous studies for RMI [START_REF] Kang | Imaging deep within a scattering medium using collective accumulation of single-scattered waves[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF][START_REF] Kwon | Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin[END_REF][START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] , except that we now have at our disposal a supplementary degree of freedom: The parameter z that controls the axial position of the focusing plane. A raw confocal image I C can be built by considering the diagonal elements of R ρρ (ρ in = ρ out ):

I C (r t , z) = R(ρ, ρ, z, t). (5)
Figure 3c shows the en-face image obtained at a given time-of-flight t for different values of z. Qualitatively, we see that the image quality strongly depends on the relative position between the coherence volume and the focusing plane. Here the presence of a highly reflecting structure, a corneal nerve, allows us to determine the parameter z that allows to match the focusing plane with the coherence volume.

Self-portrait of the focusing process.

A more quantitative and robust observable is provided by the off-diagonal coefficients of R ρρ (z, t) that enable to probe the focusing quality at any voxel. More precisely, this can be done by investigating the reflection point spread function (RPSF) defined as follows:

RP SF (∆ρ, ρ, z, t) = |R(ρ -∆ρ/2, ρ + ∆ρ/2, z, t)| 2 , (6)
This quantity derived from the off-diagonal coefficients of R ρρ , quantifies the focusing quality for each point r t = (ρ, z t ). For a medium of random reflectivity and under a local isoplanatic assumption, its ensemble average actually scales as [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF] (Supplementary Section S5):

⟨RP SF (∆ρ, ρ, z, t)⟩ ∝ |h in | 2 ∆ρ ⊛ |h out | 2 (∆ρ, ρ, z, t) (7)
where the symbols ⟨• • • ⟩ and ⊛ stand for ensemble average and convolution product, respectively. h in/out (∆ρ, ρ, z, t) is the spatial distribution of the input/output PSF along the de-scanned coordinate ∆ρ in the coherence plane at z t when trying to focus at point (ρ, z).

The RPSF can thus provide a self-portrait of the focusing process inside the cornea. Figure 3a shows the evolution of the laterally-averaged RPSF for a given time t as a function of the parameter z in the Fresnel propagator (Methods).

As expected, the focusing plane and coherence volume coincide when the RPSF extension is minimized (Fig. 3b), i.e for a defocus distance ∆z = z -z t = 0 (Fig. 3a 2 ). The estimated defocus is roughly constant over the whole thickness of the cornea. This proves that the effective index of the cornea is actually very close to the refractive index n 0 used in our propagation model (see Supplementary Section S4).

Figures 2b 1 -b 3 displays longitudinal and transverse cross-sections of the confocal image obtained after tuning the coherence volume and focusing plane at any depth (see also Supplementary Movies 1 and 2 28 ). The resolution and contrast are much better than the incoherent compound image (Supplementary Fig. S2).

In particular, the axial resolution δz c of the digital confocal image benefits from the virtual time gating and confocal filters:

δz c = 1/[2(1/δz t + 1/δz f )] ∼ 0.5 µm.
However, the image quality remains perfectible. Indeed, the RPSF still spreads well beyond the theoretical resolution cell (∼ 1 pixel) in Fig. 3a 2 . These residual aberrations originate from the lateral fluctuations of the optical index n(ρ, z) in the cornea. To demonstrate this last assertion, the transverse evolution of the focusing process can be investigated by a local assessment of the focusing quality (see Methods). A map of local RPSFs is displayed in Fig. 4a. Although the digital autofocus process provides a correct focusing quality over the whole thickness of the cornea on average, the local RPSFs exhibit important fluctuations acrossthe field-of-view. This observation is a manifestation of the 3D distribution of the optical index n(r) inside the cornea. This anisoplanic feature requires a local compensation of aberrations as we will see below.

Local Compensation of Wave Distortions.

By considering the set of autofocused reflection matrices, R ρρ (z t ), a local compensation of transverse aberrations can be performed at each depth z t [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] . It basically consists in a local analysis of wave distortions on overlapping spatial windows of size L = 18.6 µm.By exploiting a shift-shift memory effect characteristic of anisotropic scattering in the cornea [START_REF] Judkewitz | Translation correlations in anisotropically scattering media[END_REF] , one can estimate the input and output aberration phase matrices, Φ in/out (z t ) = [ϕ in/out (u in/out , ρ, z t )], between the pupil plane (u in/out ) and the medium voxels (ρ, z t ) (Methods) [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] . The result is displayed in Fig. 4b at z t = 140 µm. Strikingly, the estimated aberration laws exhibit strong phase fluctuations and vary quickly between neighboring windows. This complex feature has two origins: (i) the lateral fluctuations exhibited by the optical index inside the cornea; (i ) the imperfections of the imaging system. The latter component accounts for the difference observed between the input and output aberration transmittances (Supplementary Section S4). In fact, the input aberration phase law accumulates not only the input aberrations of the sample-arm but also those of the reference arm. The sample-induced aberrations can be investigated independently from the imperfections of the experimental set up by considering the output aberration phase matrix Φ out . The aberration phase is mainly a defocus that varies across the field-of-view due to lateral variations of the optical index.

Local shifts of the pupil function are also observed and result from a local curvature of the coherence volume with respect to the focusing plane.

The extracted aberration phase laws can be used to build transmission matrices G in/out containing the estimated impulse responses between the image voxels (ρ in/out ) and each focal plane located at depth z t (Methods). The focused R-matrix is then de-convolved by applying the phase conjugate of the transmission matrices at its input and output (Fig. 2a 3 ), such that:

R (c) ρρ (z t ) = G † out (z t ) × R ρρ (z t ) × G * in (z t ) (8)
The final image of the sample can be obtained by considering the diagonal elements of the corrected matrix R (c) ρρ : the confocal peak exhibits a spreading well beyond the diffraction limit and a background at depth due to forward multiple scattering events (Fig. 2e 1 ,f 1 ), RMI compensates for these two issues and leads to an almost ideal RPSF (Fig. 2e 2 ,f 2 ).

I M (r t ) = R (c) (ρ, ρ, z t ). ( 9
) Figures 2c 3 -
The map of final RPSFs displayed by Fig. 4c 3 shows the high focusing quality provided by RMI over the whole field-of-view at the considered depth z t = 140 µm.

The obtained three-dimensional image highlights several crucial features of the cornea: its lamellar structure induced by the collagen fibrils (Fig. 2b 3 ); (ii ) the complex network of nerves that covers the cornea; (iii ) characteristic structures of the cornea such as keratocytes and; (iv ) stromal striae whose presence is an indicator of keratoconus [START_REF] Grieve | Stromal striae: a new insight into corneal physiology and mechanics[END_REF] . Such a high-resolution image can thus be of particular 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. In that perspective, the ability of RMI in overcoming high-order aberrations and multiple scattering constitutes a paradigm shift for deep optical microscopy.

Discussion

In contrast with previous works that considered the reflection matrix at a single frequency [START_REF] Popoff | Exploiting the Time-Reversal Operator for Adaptive Optics, Selective Focusing, and Scattering Pattern Analysis[END_REF] or time-of-flight [START_REF] Kang | Imaging deep within a scattering medium using collective accumulation of single-scattered waves[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF] , the measurement of a polychromatic reflection the aberration level and scales as the number of resolution cells covered by the RPSF [START_REF] Bureau | Three-dimensional ultrasound matrix imaging[END_REF] . Compared with time-domain and scanning OCT, spectral measurement and spatial multiplexing of the wave-field provides a much better signal-to-noise ratio [START_REF] Mertz | Introduction to Optical Microscopy[END_REF] . Moreover, while a time-gated reflection matrix only allows a transverse compensation of aberrations, the polychromatic reflection matrix gives access to temporal degrees of freedom that can be exploited for compensating the axial distortions of the coherence volume. Eventually, it can be exploited for overcoming the multiple scattering limit in optical microscopy since it provides the opportunity of tailoring complex spatio-temporal focusing laws 21 required to focus light in depth.

To do so, the mapping of the refractive index will also be an important step to build accurate focusing laws inside the medium [START_REF] Chen | Multi-layer born multiple-scattering model for 3d phase microscopy[END_REF] . 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.

In that perspective, an issue we have not considered yet is medium motion during the acquisition of the reflection matrix. Of course, the assumption of a static medium is everything but true especially for in-vivo applications [START_REF] Jang | Relation between speckle decorrelation and optical phase conjugation (OPC)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skin[END_REF] . To cope with the dynamic features of the medium, two strategies can be followed. The first one is to limit the measurement time of the R-matrix at its minimum, as allowed by our device using a few illuminations. The second one is to develop algorithms that consider medium motion during the measurement of R 37 . Interestingly, temporal fluctuations of the medium's reflectivity and refractive index can provide a key information for probing the multi-cellular dynamics in optical microscopy [START_REF] Apelian | Dynamic full field optical coherence tomography: subcellular metabolic contrast revealed in tissues by interferometric signals temporal analysis[END_REF][START_REF] Scholler | Dynamic full-field optical coherence tomography: 3d live-imaging of retinal organoids[END_REF] .

Just as the concept of plane-wave imaging [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] revolutionized the field of ultrasound [START_REF] Macé | Functional ultrasound imaging of the brain[END_REF][START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF] by providing an unprecedented frame rate, our device will constitute in a near future an ideal tool for probing the 3D dynamics of tissues at a much smaller scale.

Methods

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

Sample preparation.

In the presented experiment, the corneal sample was fixed with paraformaldehyde (4% concentration).

Sampling of input and output wave-fields.

The dimension of the input pupil is D in × D in = 9 × 9 mm; the spatial sampling of input wave-fields is δ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 (δs out = 12 µm), the output wave-field is sampled at a resolution close to λ/(4N A):

δρ out = 290 nm.

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 Go) 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.

On-axis holography.

For each input wave-field, the interferogram recorded by the camera can be expressed as follows:

I(s out , u in , ω) = |E(s out , u in ,ω) + E ref (s out , u in , ω)| 2 (10)
with E and E 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(s out , u in , t) = E(s out , u in , t) t ⊛ E * (s out , u in , -t) (11)
+ E ref (s out , u in , -t) t ⊛ E * ref (s out , u in , -t) (12)
+ E ref (s out , u in , t) t ⊛ E * (s out , u in , -t) (13)
+ E(s out , u in , t) t ⊛ E * ref (s out , u in , -t) (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(s out , u 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(s out , u in , ω) = E(s out , u in , ω)E * ref (s out , u in , ω). (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 ref (s out , u in , ω) = exp i 2π λf u in .s out .
The multi-spectral reflection matrix is thus extracted using the following relation:

R(s out , u in , ω) = D(s out , u in , ω) exp -i 2π λf u in .s out . (16)

Fresnel operators

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

For this purpose, the multi-spectral reflection matrix should be first projected in the output pupil plane (u out ) by a simple 2D spatial Fourier transform:

R uu (ω) = P * us (ω) × R su (ω) (17)
where P us = [P (u, s)] is the Fourier transform operator:

P (u, s, ω) = e -j ω c 0 u•s f (18)
with f , the focal length of the MOs and c 0 the vacuum light velocity. A Fresnel phase law is then applied at the output of R uu (ω) to numerically shift the focal plane, originally located in the middle of the sample (z = 0) to any depth z:

R ρu (z, ω) = [P ⊤ uρ (ω) • F u (z, ω)] × R uu (ω) (19)
where the symbol • accounts for the Hadamard (term-by-term) product. Each column vector F 0 (z, ω) is a phase mask that accounts for the propagation of each plane wave of transverse wave vector k || = ωu/(c 0 f ) over a thickness z of an homogeneous medium of refractive index n 0 :

F 0 (u, z, ω) = e -j n 0 ω c 0 -kz z O(u) (20)
with

k z = ω c 0 n 2 0 - ||u|| 2 f 2 , (21)
F sρ (z, ω) = [P ⊤ uρ (ω) • F 0 (z, ω)] × P * us (ω) (22)
The projection between the plane wave illumination basis (u in ) and the focused basis (u in ) can also performed by means of a Fresnel propagator F ρu this time defined between the pupil plane and the each focal plane identified by their depth z:

F uρ (z, ω) = P * uρ (ω) • F 0 (z, ω) (23)

Local estimation of focusing quality

To probe the local RPSF, the field-of-view is divided at each effective depth z t into regions that are defined by their central midpoint ρ p = (x p , y p ) and their spatial extension L. A local average of the back-scattered intensity can then be performed in each region:

RP SF l (∆ρ, ρ p , z, z t ) = ⟨|R(ρ + ∆ρ/2, ρ -∆ρ/2, z, t)| 2 W L (ρ -ρ p )⟩ ρ ( 24
)
where the symbol ⟨• • • ⟩ m stands for an average over the variable m in subscript.

W L (ρ -ρ p ) = 1 for |x -x p | < L/2 and |y -y p | < L/2, and zero otherwise. In this paper, a spatial window of size L = 18.6 µm has been used to smooth out fluctuations due to the sample inhomogeneous reflectivity [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF] .

Local compensation of wave-distortions

The starting point is the time-gated reflection matrix R ρρ (z t ), obtained after tuning the focusing plane and coherence volume at each echo time t. The first step is a projection of R ρρ (z t ) in the pupil plane at input via a numerical Fourier transform:

R ρu (z t ) = P ⊤ uρ × R ρρ (z t ) (25)
An input distortion matrix is then built by performing a element-wise product between R uρ (z t ) and the phase conjugate reference matrix T uρ (ω c ) that would be obtained in absence of aberrations [START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF] (Supplementary Section S6):

D ρu (z t ) = R ρu (z t ) • P † uρ ( 26
)
A local correlation matrix C in of wave distortions is then built around each point r p = (ρ p , z t ) of the field-of-view (Supplementary Section S7). Its cofficients write:

C in (u in , u ′ in , r p ) = ⟨D(ρ out , u in , z t )D(ρ out , u ′ in , z t )W L (ρ out -ρ p )⟩ ρout (27)
Iterative phase reversal (see further) is then applied to each correlation matrix C in (r p ) [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] (Supplementary Section S8). The resulting input phase laws, Φ in (z t ) =

[ϕ in (u in , ρ p , z t )], are used to compensate for the wave distortions undergone by the incident wave-fronts:

R ′ ρρ (z t ) = {R ρu (z t ) • exp [-iϕ in (z t )]} × P * uρ ( 28
)
The corrected matrix R ′ ρρ is only intermediate since phase distortions undergone by the reflected wave-fronts remain to be corrected.

To that aim, R ′ ρρ (z t ) is now projected in the pupil plane at output:

R ′ uρ (z t ) = P uρ × R ′ ρρ (z t ). ( 29
)
An output distortion matrix is then built:

D uρ (z t ) = P * uρ • R ′ uρ (z t ) (30)
From D uρ , one can build a correlation matrix C out for each point r p :

C out (u out , u ′ out , r p ) = ⟨D(u out , ρ in , z t )D * out (u ′ out , ρ in , z t )W L (ρ in -ρ p )⟩ ρ in (31)
The IPR algorithm described further is then applied to each matrix C out (r p ).

The resulting output phase laws, Φ out (z t ) = [ϕ out (u out , ρ p , z t )], are leveraged to compensate for the residual wave distortions undergone by the reflected wavefronts:

R (c) ρρ (z t ) = P † uρ × exp [-iΦ out (z t )] • R ′ uρ (z t ) ( 32
)
Iterative phase reversal algorithm.

The IPR algorithm is a computational process that provides an estimator of the phase of the transmittance that links each point u of the pupil plane with each voxel r p = (ρ p , z t ) of the cornea volume [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] . To that aim, the correlation matrix C computed over the spatial window W L centered around a given point (ρ p , z t ) is considered (Eqs. 27 and 31). Mathematically, the algorithm is based on the following recursive relation:

Φ (n) in/out (ρ p , z t ) = arg C(ρ p , z t ) × exp iΦ (n-1) in/out (ρ p , z t ) ( 33
)
where

Φ (n)
in/out is the estimator of the transmittance phase at the n th iteration of the phase reversal process. Φ (0) in/out is an arbitrary wave-front that initiates the process (typically a flat phase law) and ϕ in/out = lim n→∞ ϕ (n) in/out is the result of IPR.

Transmission matrices.

The transmission matrices G in/out used to deconvolve the focused R-matrix (Eq. 8) can be deduced from the estimated aberration phase laws Φ in/out as follows (Eqs. 28 and 32):

G in/out (z t ) = P uρ • exp iΦ in/out (z t ) × P † uρ . (34

S3. MULTI-SPECTRAL REFLECTION MATRIX

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

The wave field E s (s out , u in , ω) reflected by the sample arm in the camera plane can be expressed as follows:

E s (s out , u in , ω) = S(ω) G(s out , ρ s , z s , ω)γ(ρ s , z s )E in (u in , ρ s , z s , ω)dρ s dz s . (S1)
S(ω) is the amplitude of light source at frequency ω. G(s out , ρ s , z s ), the Green's function between the sample mapped by the vector (ρ s , z s ) and the CCD sensors identified by s out . γ(ρ s , z s ) represents the sample reflectivity. E in (u in , ρ s , z, ω)

describes the incident wave-field that can be expressed as follows:

E in (u in , ρ s , z s , ω) = S(ω)F(u in , z s , ω)T (u in , ρ s , z) exp -i 2π λf u in .ρ s (S2)
where F(u in , z s , ω) is the Fresnel phase law that describes plane wave propagation inside an homogeneous medium of effective index n (Eq. 20 in the accompanying paper). The transmission matrix T = [T (u in , ρ s , z)] accounts for the wave distortions induced by the fluctuations of the optical index inside the medium.

In the reference arm, a mirror is placed in the focal plane of the MO and displays a uniform reflectivity: γ(ρ m , z m ) = γ m δ(z m ), with δ the Dirac distribution and γ m the mirror surface reflectivity. The reference wave-field is thus given by:

E ref (s out , u in , ω) = γ m S(ω) G ref (s out , ρ m )T ref (u in , ρ m ) exp -i 2π λf u in .ρ m dρ m . (
T ref (u in ) = dρH ref (ρ) exp -i 2π λf u in • ρ . (S4)
Under the isoplanatic assumption, Equation S3 thus simplifies into:

E ref (s out , u in , ω) = γ m S(ω)T ref (-u in )T ref (u in ) exp i 2π λf u in .s out (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 ref (s out , u in , ω) = γ m S(ω) exp i 2π λf u in .s out . (S6)
In the following, we will not make this assumption and will consider the more general expression of E ref given in Eq. S5.

The coefficients of the multi-spectral matrix R su (ω) are recorded by isolating the interference between the sample beam, E s , and the reference beam, E ref (Eqs. 15 and 16 of the accompanying paper):

R(s out , u in , ω) = E(s out , u in , ω)E * ref (s out , u in , ω) exp i 2π λf u in .s out (S7)
Using Eqs. S1, S2, S5, the last equation can be rewritten as follows:

R(s out , u in , ω) = γ m |S(ω)| 2 G(s out , ρ s , z s , ω)γ(ρ s , z s ) × T (u in , ρ s , z s )F(u in , z s , ω) × T * ref (-u in )T * ref (u in ) exp -i 2π λf u in .ρ s dρ s dz s . (S8)
One can already notice from this expression that the aberrations undergone by the reference wave-field (T * ref (-u in )T * ref (u in )) emerge at the input of the recorded reflection matrix.

S4. FOCUSED REFLECTION MATRIX

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

First, a spatial Fourier transform over the camera pixels s out leads to the reflection matrix R uu in the pupil basis (Eqs. 17 and 18 of the accompanying paper):

R(u out , u in , ω) = γ m |S(ω)| 2 T (u out , ρ s , z)F(u out , z s , ω)γ(ρ s , z s ) × T (u in , ρ s , z)F(u in , z s , ω) × T * ref (-u in )T * ref (u in ) (S9) × exp -i 2π λf (u in + u out ).ρ s dρ s dz s . (S10)
As for incident light (Eq. S2), the return path is decomposed in the plane wave basis as the product between a Fresnel phase law F(u out , z s , ω), accounting for free-space wave propagation in an homogeneous medium of refractive index n, and the transfer function T (u in , ρ s , z s ) that grasps the wave distortions induced by the refractive index fluctuations such that:

T (u out , ρ s , z s )F(u out , z s , ω) = G(s out , ρ s , z s ) exp -i 2π λf u out .s out dρ s . (S11)
Numerical focusing at depth z (Eqs. 1 and 4 of the accompanying paper) then consists in compensating wave diffraction by applying the phase conjugate of the Fresnel propagator for refractive index n 0 at input and output before a spectral Fourier transform:

R(ρ out , ρ in , z, t) = γ m u in uout dω|S(ω)| 2 exp (iωt) F * 0 (u out , z, ω)R(u out , u in , ω)F * 0 (u in , z, ω) × exp i 2π λf (u in • ρ in + u out • ρ out ).ρ s
Injecting Eq. S9 leads to the following expression for the coefficients of R ρρ (z, ω):

R(ρ out , ρ in , z, t) =γ m u in uout dω|S(ω)| 2 exp (iωt) × T (u out , ρ s , z)γ(ρ s , z s , ω)T (u in , ρ s , z s )T * ref (-u in )T * ref (u in ) × F(u out , z s , ω)F * 0 (u out , z, ω)F(u in , z s , ω)F * 0 (u in , z, ω) × exp -i 2π λf [u in .(ρ s -ρ in ) + u out .(ρ s -ρ out )] dρ s dz s . (S12)
The positions of the coherence volume and focusing plane are determined by the cancellation of the Fresnel phase laws,

F(u in , z, ω)F * 0 (u in , z s , ω)F(u out , z, ω)F * 0 (u out , z s , ω) = exp i ω c 0 n 2 0 - ||u in || 2 f 2 + n 2 0 - ||u out || 2 f 2 -2n 0 z × exp -i ω c 0 n 2 - ||u in || 2 f 2 + n 2 - ||u out || 2 f 2 z s . (S13)
Under the paraxial approximation, these Fresnel phase laws can be developed as follows:

F(u in , z, ω)F * 0 (u in , z s , ω)F(u out , z, ω)F * 0 (u out , z s , ω) = exp -2i nω c 0 z s exp -i ω c 0 ||u in || 2 2f 2 + ||u out || 2 2f 2 z n 0 - z s n , (S14)
The cancellation of the first phase term defines the real position of the coherence volume z s = z t = c 0 t/(2n) that appears at an effective depth z 0 = c 0 t/(2n 0 ) = (n/n 0 )z t (Fig. S3a). Previous expression of R(ρ out , ρ in , z, t) (Eq. S12) can be rewritten as follows:

R(ρ out , ρ in , z, t) =γ m f (t -2nz t /c) t ⊛ (S15) ω u in uout T (u out , ρ s , z)γ(ρ s , z s ) (S16) × T (u in , ρ s , z s )T * ref (-u in )T * ref (u in ) exp (iωt) × exp -i ω c 0 ||u in || 2 2f 2 + ||u out || 2 2f 2 z n 0 - z s n × exp -i 2π λf [u in .(ρ s -ρ in ) + u out .(ρ s -ρ out )] dρ s dz s . (S17)
with f (t) = dω|S(ω)| 2 e iωt , the time response of the microscope. The symbol ⊛ stands for convolution over variable t. The position z f of the focusing plane is obtained when the parabolic phase term cancels in the previous expression, that is to say for z s = z f = (n 0 /n)z t . The apparent defocus induced by the mismatch between n and n 0 is thus equal to:

∆z = z f -z 0 = z t (n 0 /n -n/n 0 ) (S18)
An index mismatch thus implies a defocus distance ∆z that increases linearly with z t (Fig. S3a). In the present study, the estimated defocus is roughly constant with 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 ∆z can be applied in post-processing in order to make coincide the coherence and focusing planes. 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 ρρ (z t ) can be derived as follows:

R(ρ out , ρ in , z t ) =γ m u in uout T out (u out , ρ s , z t )γ(ρ s , z t )T in (u in , ρ s , z t ) × exp i 2π λf [u in .(ρ s -ρ in ) + u out .(ρ s -ρ out )] dρ s . (S19) with T out (u out , ρ s , z t ) ≡ T (u out , ρ s , z t )
and

T in (u in , ρ s , z t ) ≡ T (u in , ρ s , z t )T * ref (-u in )T * ref (u in ).
While the output transmission matrix T out corresponds to the sample transmission matrix T, the input transmission matrix T in grasps both the sample and reference arm aberrations.

The Fourier transform of the transmission coefficients T in/out in Eq. S19 provide local PSFs, H in/out , such that:

H in/out (ρ, ρ s , z t ) = u T in/out (u, ρ s , z t ) exp i 2π λf u.ρ . (S20)
The time-gated reflection matrix can be rewritten as follows:

R(ρ out , ρ in , z t ) = dρ s H out (ρ s -ρ out , ρ s , z t )γ(ρ s , z)H in (ρ s -ρ in , ρ s , z t ) (S21)
The latter expression can be recast as a function the impulse responses G in/out (ρ s , ρ in/out , z t ) between input/output focusing points and points ρ s mapping the sample. Both quantities are actually linked as follows:

G in/out (ρ s , ρ in/out , z t ) = H in/out (ρ s -ρ in/out , ρ s , z t ) (S22)
Injecting the last expression into Eq. S21 leads to the following expression:

R(ρ out , ρ in , z t ) = G out (ρ out , ρ s , z t )γ(ρ s , z t )G in (ρ in , ρ s , z t ) (S23)
Under a matrix formalism, the last expression can be rewritten as follows:

R ρρ (z t ) = G out (z t ) × Γ(z t ) × G ⊤ in (z t ), (S24)
Γ describes the scattering process inside the sample. Under a single scattering assumption, this matrix is diagonal. Its coefficients then correspond the sample reflectivity γ(ρ s , z t ). G in and G out are the input and output focusing matrices.

Their coefficients, G in/out (ρ in/out , ρ s , z t ), describe the transverse amplitude distri-bution of the focal spot when trying to focus at point (ρ in/out , z t ).

S5. REFLECTION POINT SPREAD FUNCTION

As mentioned in the accompanying paper, the off-diagonal cofficients of R ρρ (z, t)

enable to probe the focusing quality at any voxel by investigating the reflection point spread function (RPSF, Eq. 6). To express theoretically the latter quantity, a local isoplanatic assumption shall be made. This hypothesis implies that the PSFs H in/out (ρ s -ρ in/out , ρ s , z, t) = h in/out (ρ s -ρ in/out , ρ p , z p , t).

(S25)

The second assumption is to consider the medium reflectivity γ(ρ s , z) as random:

⟨γ (ρ 1 , z) γ * (ρ 2 , z)⟩ = |γ| 2 δ (ρ 2 -ρ 1 ) , (S26)
By combining those assumptions with Eq. S21, the ensemble average of RP SF (∆ρ, ρ p , z, t) (Eq. 7) can be expressed as follows:

⟨RP The physical meaning of this last equation is the following: Each distorted wavefield corresponds to the diffraction of a virtual source synthesized inside the medium modulated by the transmittance T out of the sample between the focal and pupil planes. Each virtual source is spatially incoherent due to the random reflectivity of the medium, and its size is governed by the spatial extension of the input focal spot. The idea is now to smartly combine each virtual source to generate a coherent guide star and estimate T out independently from the sample reflectivity.

incomplete convergence towards the covariance matrix leads to a phase error δϕ made on our estimation of each aberration phase law. The variance of this error scales as follows:

|δϕ out (u)| 2 ∼ 1 N L C in (S42)
with C in , the coherence factor that is a direct indicator of the focusing quality [START_REF] Mallart | Adaptive focusing in scattering media through soundspeed inhomogeneities: The van Cittert Zernike approach and focusing criterion[END_REF] , ranging from 0 for a fully blurred guide star to 2/9 for a diffraction-limited focal spot [START_REF] Silverstein | Ultrasound scattering model: 2-d cross-correlation and focusing criteria-theory, simulations, and experiments[END_REF] . On the one hand, this scaling of the phase error with N W explains why each spatial window should be large enough to encompass a sufficient number of independent realizations of disorder. On the other hand, it should be small enough to grasp the spatial variations of aberration phasr laws across the field-of-view. A compromise has thus to be found between these two opposite requirements. It has led us to to choose spatial windows of size L = 18.6 µm to ensure the convergence of the IPR process. Note that, contrary to a recent study performed in an extremely opaque cornea [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] , a multi-scale approach of wave distortions is not necessary here since the cornea under study displays smoother fluctuations of optical index over larger characteristic length scales. As previous works dealing with the CLASS algorithm [START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF] , 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.

[1] V. Ntziachristos, Going deeper than microscopy: The optical imaging frontier in biology, Nat. Methods 7, 603 (2010).

[2] S. Gigan, Optical microscopy aims deep, Nat. Photonics 11, 14 (2017).

FIG. 1 .

 1 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 (u in ) in each arm (a 1 ). The sample placed in the focal plane of the first MO (MO 1 , NA=0.8) is thus illuminated by a set of plane waves at each frequency of the light source bandwidth (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 (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 2 , NA=0.8). b, At each frequency ω, for each input wave-front u in , the interferogram I(s out ) (b 1 ) recorded by each pixel s out of the camera provides one column of the spectral reflection matrix R su (ω) = [R(s out , u in , ω)] (b 2 ).

FIG. 2 .

 2 FIG. 2. From Holoscopy to Matrix Imaging. a, Imaging Methods. a 1 , 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 r out mapping the sample (red). a 2 , 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) at each depth z t . a 3 , 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 = 150 µm and 275 µ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 3 , blue and red arrows design some corneal nerves and keratocytes , respectively; in panel d 3 , stroma striae are highlighted by green dashed lines. e, RPSF at depth z = 150 µm [scale bar: 3 µm] for DCM (e 1 -f 1 ) and RMI (e 2 -f 2 ) images. The radial evolutions of these RPSFs are compared in panel (e 3 ) [DCM: blue; RMI: orange]. f, Same as in panel (e) but at depth z = 275 µm.

FIG. 3 .

 3 FIG. 3. Auto-focusing process guided by the reflection point spread function. a Evolution of the RPSF versus the defocus distance ∆z for a fixed coherence volume. The transverse distribution of the RPSF is shown for several values of defocus (a 1 )-(a 3 ). The evolution of its radial average is displayed in panel (a 4 ). b. Relative position of the focusing plane (dash-dotted line) and coherence volume (red layer) for the different values of defocus ∆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 ∆z = -4, 0, and +4 µm. The considered coherence volume is located at the effective depth z t = 140 µm in the cornea.

FIG. 4 .

 4 FIG. 4. Local Adaptive Focusing by Matrix Analysis of Wave Distortions. a, Transverse map of local RPSFs after the auto-focusing process (a 1 ) and after the aberration matrix compensation of wave distortion (a 2 ) at z t = 140 µm Each RPSF is displayed over a de-scan area of 6 × 6 µm 2 . b, Transverse map of input (b 1 ) and output (b 2 ) aberration phase laws, ϕ in/out , estimated from the pupil plane at the same depth.

  Figures 2c 3 -e 3 display the corresponding longitudinal and transverse cross-sections obtained by RMI (see also Supplementary Movies 1, 2 and 3 28 ). The comparison with the confocal image [Figs. 2c 2 -e 2 ] shows a clear gain in contrast. The resolution improvement can be assessed by examining the RPSF. While, at the previous step,

  matrix 22 allowed us to realize in post-processing: (i) a 3D confocal image of the sample reflectivity on millimetric volumes (0.1 mm 3 = 10 9 voxels) in an ultrafast acquisition time (1 s) ; (ii) a local compensation of aberrations which usually prevent deep imaging. The required number of input wave-fronts depends on

  the longitudinal component of the wave vector, and O(u), the finite pupil support: O(u) = 1 for ||u|| < f N A and zero elsewhere. Each reflection matrix R ρu (z, ω) = [R(ρ out , u in , z, ω)] connects each output virtual focusing point r out = (ρ out , z) to each input illumination u in at frequency ω. The combination of Eqs. 17, 19 and 20 leads us to define the Fresnel operator F sρ (z, ω) = [F (s out , ρ out , z, ω)] that enables the projection of the optical data from the camera sensors (s out ) to any focal plane (ρ out ) of expected depth z (Eq. 1):

  FIG. S1. Detailed experimental set up. P: Polarizers; A: Analyzer; QW: quarter wave plates.

FIG. S2 .

 S2 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).

  S3) with G ref (s out , ρ m , ω), the Green's functions between the focal plane of the MO (ρ m ) and the CCD sensors (s out ) and T ref (u in , ρ 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 ref can be replaced by a spatiallyinvariant impulse response H ref between the focal plane and the CCD sensors: G ref (s out , ρ m ) = H ref (ρ m +s out ). Under the same hypothesis, the transfer function T ref (u in , ρ m ) becomes an aberration transmittance T ref (u in ), defined as the Fourier transform of the reference arm point spread function H ref :

  z t . It thus means that the cornea displays an effective optical index n ∼ n 0 = 1.33 and that the observed defocus rather originates from the different lengths between the sample and reference arms.

H

  in/out are locally invariant by translation. This leads us to define local spatiallyinvariant PSFs h in/out around each central midpoint ρ p at each time-of-flight such that:

  SF (∆ρ, ρ p , z, t)⟩ = |γ| 2 × |h in | 2 ∆ρ ⊛ |h out | 2 (∆ρ, ρ p , z p , t).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. The output distortion matrix can be built by first projecting the time-gated are extracted and stored in a local distortion matrix D ′ out (ρ p , z t ):D ′ (u out , ρ in , ρ p , z t ) = D(u out , ρ in , z t ) W L (ρ in -ρ p ), (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 out (u out , ρ p , z t ) around each point (ρ p , z t ), such thatT out (u out , ρ s , z t ) ≃ T out (u out , ρ p , z t). This transmittance is the Fourier transform of the local PSF h out (ρ, ρ p , z t ):T out (u out , ρ p , z t ) = h out (ρ, ρ p , z) exp 2π λf u out • ρ dρ (S34)Under this assumption, Eq.S31 can be rewritten as follows:D ′ (u out , ρ in , ρ p , z t ) = T out (u out , ρ p , z t )

	(S35)
	transmittance
	(S27)
	S6. LOCAL DISTORTION MATRIX

× γ(ρ + ρ out , z t )h in (ρ, ρ p , z t ) exp -i 2π λf u in • ρ dρ virtual source

.

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

Funding Information. The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement nos. 610110 and 819261, HELMHOLTZ* and REMINISCENCE projects, respectively). This project has

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

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

Supplementary Information

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; An analyzer A1 allows us to project the sample and reference beams on the same polarization and make them interfere in the camera plane.

focused reflection matrix R ρρ (z t ) in the pupil plane at output:

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:

or in terms of matrix coefficients,

Injecting Eqs. S20 and S21 into the last equation yields

In previous papers [START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] , we showed that the distortion matrix D highlights spatial correlations of the reflected wave-field induced by the shift-shift memory effect [START_REF] Judkewitz | Translation correlations in anisotropically scattering media[END_REF][START_REF] Osnabrugge | Generalized optical memory effect[END_REF] :

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-ofillumination should be subdivided into overlapping regions [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] that are defined by their central midpoint (ρ p , z t ) and their lateral extension L. All of the distorted components associated with focusing points ρ in located within each region

S7. CORRELATION OF WAVE DISTORTIONS

To do so, the correlation matrix C out = D out D † out is an excellent tool. Its coefficients write as follows

The matrix C out (ρ p , z t ) can be decomposed as the sum of its ensemble average, the covariance matrix ⟨C out ⟩ (ρ p , z t ), and a perturbation term δC out (ρ p , z t ):

The intensity of the perturbation term scales as the inverse of the number N L = (L/δρ 0 ) 2 of resolution cells in each sub-region [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF][START_REF] Lambert | Ultrasound matrix imaging -part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] :

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

Under assumptions of local isoplanicity (Eqs. S25 and S35) and random reflectivity,

The coefficients of the covariance matrix can be expressed as follows [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] :

or in terms of matrix coefficients,

where the symbol * stands for correlation product over variable u. C H is a reference correlation matrix that would be measured in an homogeneous cornea for a virtual reflector whose scattering distribution corresponds to the output focal spot intensity |h in (ρ s , ρ p , z t )| 2 . The covariance matrix ⟨C in ⟩ (ρ p , z t ) thus corresponds to the same experimental situation but for a virtual reflector embedded into the heterogeneous cornea under study.

S8. ITERATIVE PHASE REVERSAL

An estimator exp[iΦ out (ρ p , z t )] of the local aberration transmittance T out (ρ p , z t ) can be extracted by applying an iterative phase reversal algorithm to C out (Eq. 34 of the accompanying paper). It mimics an iterative time reversal process on the virtual reflector described above but imposes a constant amplitude for the timereversal invariant. This iterative phase reversal (IPR) process converges towards a wavefront that maximizes the coherence of the wave-field reflected by the virtual reflector [START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] .

This IPR process assumes the convergence of the correlation matrix C out (Eq. S37) towards its ensemble average ⟨C out ⟩, the covariance matrix [START_REF] Lambert | Ultrasound matrix imaging -part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF][START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] . In fact, this convergence is never fully realized and C out should be decomposed as the sum of this covariance matrix ⟨C out ⟩ and the perturbation term δC out (Eq. S37). This