Ulysse Najar

Victor Barolle

Paul Balondrade

Mathias Fink

A Claude Boccara

Alexandre Aubry

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

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

Multiple scattering of waves concerns many domains of physics, ranging from optics or acoustics to solid-state physics, seismology, medical imaging, or telecommunications. In an inhomogeneous medium where the refractive index n depends on the spatial coordinates r, several physical parameters are relevant to characterize wave propagation: (i) the scattering mean free path ℓ s , which is the average distance between two successive scattering events;

(ii) the transport mean free path ℓ t , which is the distance after which the wave has lost the memory of its initial direction. For a penetration depth z smaller than ℓ s , ballistic light is predominant and standard focusing methods can be employed; for z > ℓ s , multiple scattering events result in a gradual randomization of the propagation direction before reaching the diffusive regime for z > ℓ t . Although it gives rise to fascinating interference phenomena such as perfect transmission 1,2 or Anderson localization 3,4 , multiple scattering still represents a major obstacle to deep imaging and focusing of light inside complex media 5,6 .

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

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

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

Inspired by previous works 25,26 , full-field optical coherence tomography (FFOCT) 27,28 will be used here to record the R-matrix. In FFOCT, the incident wave-field is temporallyand spatially-incoherent. It enables, by means of low coherence interferometry, a parallel acquisition of a time-gated confocal image 29 at a much better signal-to-noise ratio than a traditional point scanning scheme for equal measurement time and power [START_REF] Mertz | Introduction to Optical Microscopy[END_REF] . By splitting the incident wave-field into two laterally-shifted components, a de-scanned measurement of R can be performed without a tedious raster scanning of the field-of-view 20 .

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

The experimental proof-of-concept presented in this paper is performed on a human ex-vivo cornea that we chose deliberately to be extremely opaque. Its overall thickness is of 10ℓ s . FFOCT shows an imaging depth limit of 2ℓ s due to aberration and scattering.

Strikingly, RMI enables to recover a full 3D image of the cornea at a resolution close to λ/4

(∼ 230 nm) and a penetration depth enhanced by, at least, a factor five.

Measuring a De-scan Reflection Matrix

Our approach is based on a de-scanned measurement of the time-gated reflection matrix R from the scattering sample. Inspired by time-domain FFOCT 27,28 , the corresponding set up is displayed in Fig. 1a. It consists in a Michelson interferometer with microscope objectives in both arms (Fig. 1a). In the first arm, a reference mirror is placed in the focal plane of a microscope objective (MO). The second arm contains the scattering sample to be imaged. Because of the broad spectrum of the incident light, interferences occur between the two arms provided that the optical path difference through the interferometer is close to zero.

The length of the reference arm determines the slice of the sample (coherence volume) to be imaged and is adjusted in order to match with the focal plane of the MO in the sample arm.

The backscattered light from each voxel of the coherence volume can only interfere with the light coming from the conjugated point of a reference mirror. The spatial incoherence of the light source actually acts as a physical confocal pinhole (Fig. 1c). All these interference e. Corresponding RPSFs after the matrix imaging process. f. Longitudinal (x,z) section of the volumetric image at the end of the matrix imaging process. g. Schematic of a healthy human cornea. Each image is normalized at each depth by its averaged intensity.

To overcome the multiple scattering phenomenon, one should go beyond a simple confocal image and record the cross-talk between the camera pixels. Experimentally, it consists in measuring the reflection matrix R associated with the sample (Fig. 3a). Interestingly, this can be done by slightly modifying the illumination scheme of the FFOCT device, as displayed in Fig. 1a. The incident wave-fields are still identical in each arm but are laterally shifted with respect to each other by a transverse position ∆ρ in . Their spatial incoherence now acts as a de-scanned pinhole that gives access to the cross-talk between distinct focusing points operator associated with a virtual guide star that results from a coherent average of all the descanned focal spots displayed in b (Supplementary Section S7). f. IPR is then applied (Methods).

The resulting wave-front compensates for aberrations and scattering inside the medium to produce a sharper guide star. It provides an estimation of one column of T in corresponding to the common mid-point r p of the input focal spots considered in panel a.

(Fig. 1b). The interferogram recorded by the camera (Fig. 1d) directly provides one line of the reflection matrix R in de-scanned at input (Figs. 3b andc), such that

R in (∆ρ in , ρ out , z) = R(ρ out + ∆ρ in , ρ out , z), (1)
with R = [R(ρ in , ρ out , z)], the reflection matrix expressed in the canonical basis. Its coefficients R(ρ in , ρ out , z) correspond to the response of the medium at depth z between points ρ in and ρ out in the source and camera planes (Fig. 3a). Scanning the relative position ∆ρ in is equivalent to recording the canonical R-matrix diagonal-by-diagonal (see Figs. 3a andc).

However, while a raster scan (column-by-column acquisition) of R requires to illuminate the sample over a field-of-view Ω with N = (Ω/δ 0 ) 2 input wave-fronts 20,[START_REF] Badon | [END_REF]32 , the de-scanned basis allows a much smaller number of field measurements.

This sparsity can be understood by expressing theoretically the de-scan matrix R in (Supplementary Section S3):

R in (∆ρ in , ρ out , z) = Ω dρ s H in (ρ s + ∆ρ in , ρ in , z)γ(ρ s + ρ out , z)H out (ρ s , ρ out , z) ( 2
)
where γ is the sample reflectivity. The de-scanned elements allow us to go far beyond standard confocal imaging. In particular, they will be exploited to unscramble the local input and output PSFs in the vicinity of each focal point. As a preliminary step, they can also be used to quantify the level of aberrations and multiple scattering. In average, the de-scanned intensity,

I(∆ρ in , ρ out , z) = |R in (∆ρ, ρ out , z)| 2
, can actually be expressed as the convolution between the incoherent input and output PSFs 33 :

⟨I(∆ρ in , ρ out , z)⟩ ∝ |H in | 2 ∆ρ in ⊛ |H out | 2 (∆ρ in , ρ out , z) (3)
where the symbol ⊛ stands for correlation product and ⟨• • • ⟩ for ensemble average. This quantity will be referred to as RPSF in the following (acronym for reflection PSF). Figure 1e displays examples of RSPF extracted in depth of the opaque cornea. The spatial extension δ R of the RPSF indicates the focusing quality and dictates the number M of central lines of R in (z) that contain the relevant information for imaging:

M ∼ (δ R /δ 0 ) 2 (4)
with δ 0 ∼ λ/(4N A), the confocal maximal resolution of the imaging system. For a field-ofview much larger than the spatial extension of the RPSF (Ω ≫ δ R ), the de-scanned basis is thus particularly relevant for the acquisition of R (M ≪ N ).

Quantifying the Focusing Quality

Figure 2d shows the depth evolution of the RPSF. It exhibits the following characteristic shape: a distorted and enlarged confocal spot on top of a diffuse background 33 . The former component is a manifestation of aberrations; the latter contribution is due to multiple scattering. Figure 2d clearly highlights two regimes. In the epithelium (z < 70 µm), the confocal component is predominant and the image of the cornea is reliable although its resolution is affected by aberrations (Fig. 2b 1 ). Beyond this depth, the multiple scattering background is predominant and drastically blurs the image (Fig. 2b 3 ). The axial evolution of the confocal-to-multiple scattering ratio enables the measurement of the scattering mean free path ℓ s 34 (Supplementary Section S2). We find ℓ s ∼ 35 µm in the stroma (Fig. 2g), which confirms the strong opacity of the cornea. The penetration depth limit thus scales as ℓ s . This value is modest compared with theoretical predictions 8 (∼ 4ℓ s ) but is explained by the occurrence of strong aberrations at shallow depths, partially due to the index mismatch at the cornea surface (Fig. 2d).

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

Iterative Phase Reversal of Wave Distortions

To that aim, we will exploit and extend the distortion matrix concept introduced in a previous work 19 . Interestingly, a Fourier transform over the coordinate ∆ρ in of each de-scanned wave-field, R in (∆ρ in , ρ out , z), actually yields the wave distortions seen from the input pupil plane (Fig. 3d) :

D in (z) = F × R in (z) (5)
where F denotes the Fourier transform operator, F (u, ∆ρ) = exp (-i2πu.∆ρ/λf ), λ the central wavelength and f the MO focal length. D in (z) = [D(u in , ρ out, , z)] is the distortion matrix that connects any voxel (ρ out , z) in the field-of-view to wave-distortions in the input pupil plane (u in ).

As expected in most of biological tissues, this matrix exhibits local correlations that can be understood in light of the shift-shift memory effect [START_REF] Judkewitz | [END_REF]37 : Waves produced by nearby points inside an anisotropic scattering medium generate highly correlated random speckle patterns in the pupil plane. Figure 1 illustrates this fact by displaying an example of distortion matrix (Fig. 1g) and reshaped distorted wave-fields for different points (ρ out , z) (Fig. 1h). A strong similarity can be observed between distorted wave-fronts associated with neighboring points but this correlation tends to vanish when the two points are too far away.

The next step is to extract and exploit this local memory effect for imaging. To that aim, a set of correlation matrices C in (r p ) shall be considered between distorted wave-fronts in the vicinity of each point r p in the field-of-view (Methods). Under the hypothesis of local isoplanicity, each matrix C in (r p ) is analogous to a R-matrix associated with a virtual reflector synthesized from the set of output focal spots 21 (see Fig. 3e and Supplementary Section S6).

In this fictitious experimental configuration, an iterative phase-reversal (IPR) process can be performed to converge towards the incident wave front that focuses perfectly through the heterogeneities of the medium onto this virtual guide star (see Fig. 3f and Methods).

IPR repeated for each point r p yields an estimator Tin of the time-gated transmission matrix, T in = F × H in . Its digital phase conjugation enables a local compensation of aberration and forward multiple scattering. An updated de-scanned matrix can then be built:

R in = F † × T * in • D in (6)
where the symbol † stands for transpose conjugate and • for the Hadamard product. The same process can be repeated by exchanging input and output to estimate the output transmission matrix T out (Methods).

Multi-Scale Analysis of the Distortion Matrix

A critical aspect of RMI is the choice of the spatial window over which wave distortions shall be analyzed. On the one hand, the isoplanatic assumption is valid for low-order aberrations that are associated with extended isoplanatic patches. On the other hand, forward multiple scattering gives rise to high-order aberrations that exhibit a coherence length that decreases with depth until reaching the size of a speckle grain beyond ℓ t 36 .

However, each spatial window should be large enough to encompass a sufficient number of independent realizations of disorder 38 . Indeed, the bias of our T-matrix estimator scales as follows (see Supplementary Section S9):

|δT (u, r p )| 2 ∼ 1/(C 2 N W ) (7)
with N W the number of resolution cells in each spatial window. C is a coherence factor that is a direct indicator of the focusing quality 39 , ranging from 0 for a fully blurred guide star to 2/9 for a diffraction-limited focal spot 40 .

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

As discussed further, the end of the process is monitored by the memory effect that shall exhibit our T-matrix estimator (Supplementary Section S12).

Transmission Matrix and Memory Effect

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

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

The memory effect is also a powerful tool to monitor the convergence of the IPR process.

When the spatial window is too small (3×3 µm 2 ), IPR provides a spatially-incoherent estimator T and leads to a bucket-like image (Supplementary Figure S8). This observable thus indicates when the convergence towards T is fulfilled or when the algorithm shall be stopped.

Deep Volumetric Imaging

Eventually, the estimated T-matrices can be used to compensate for local aberrations over the whole field-of-view. To that aim, a digital phase conjugation is performed at input and output (Eq. 6). The comparison between the initial and resulting images (Figs. 4a,b)

demonstrates the benefit of a local compensation of aberration and scattering. The drastic gain in resolution and contrast provided by RMI enables to reveal a rich arrangement of biological structures (cells, striae, etc.) that were completely blurred by scattering in the initial image. For instance, a stromal stria, indicator of keratoconus 43 , is clearly revealed on the RMI B-scan (Fig. 2f) while it was hidden by the multiple scattering fog on the initial image (Fig. 2c). The B-scan shows that RMI provides a full image of the cornea with the recovery of its different layers throughout its thickness (350 µm ∼ 10ℓ s , see also

Supplementary Movies).

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

Discussion

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

First, the proposed IPR algorithm outperforms iterative time reversal processing 19 for local compensation of aberrations in scattering media because it can evaluate the focusing laws over a larger angular domain (Supplementary Figure S3). Second, the bias of our T-matrix estimator has been expressed analytically (Eq. S40) as a function of a coherence factor that grasps the blurring effect of aberrations and multiple scattering. This led us to define a multi-scale strategy for matrix imaging with a fine monitoring of its convergence based on the memory effect. The latter observable is a real asset as it provides an objective criterion to: (i) optimize the resolution of our T-matrix estimator (Supplementary Section S12);

(ii) compare our approach with alternative methods such as the CLASS algorithm 20,23,24 (Supplementary Section S18). Our multi-scale process enables us to target isoplanatic areas more than four times smaller than CLASS. Interestingly, those two approaches are based on the maximization of different physical quantities: the confocal intensity for CLASS; the coherence of the wave-field induced by a virtual guide star for IPR. Hence they are, in principle, perfectly complementary and could be advantageously combined in the future.

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

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

To conclude, this study is a striking illustration of a pluri-disciplinary approach in wave physics. A passive measurement of the R-matrix is indeed an original idea coming from seismology 50 . The D-matrix is inspired by stellar speckle interferometry in astronomy 51 .

The T-matrix is a concept that has emerged both from fundamental studies in condensed matter physics 52 and more applied fields such as MIMO communications 53 and ultrasound therapy 12 . The emergence of high-speed cameras and the rapid growth of computational capabilities now makes matrix imaging mature for deep in-vivo optical microscopy.

Methods

Experimental set up

The full experimental setup is displayed in Supplementary Figure S1.

Cornea

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

Experimental procedure

The experiment consists in the acquisition of the de-scanned reflection matrix R in . To that aim, an axial scan of the sample is performed over the cornea thickness (350 µm) with a sampling of 2 µm (i.e 185 axial positions). For each depth, a transverse scan of the de-scanned position ∆ρ in is performed over a 2.9 × 2.9 µm 2 area with a spatial sampling δ 0 = 230 nm (that is to say 169 input wave-fronts instead of 10 6 input wave-fronts in a canonical basis).

For each scan position (∆ρ, z), a complex-reflected wave field is extracted by phase shifting interferometry from four intensity measurements. This measured field is averaged over 5 successive realisations (for denoising). The integration time of the camera is set to 5 ms.

Each wave-field is stored in the de-scanned reflection matrix R in = [R in (∆ρ in , ρ out )] (Fig. 1).

The duration time for the recording of R in is of ∼ 30 s at each depth. The post-processing of the reflection matrix (IPR and multi-scale analysis) to get the final image took only a few minutes on Matlab. The experimental results displayed in Fig. 4 and 5 at a single depth z = 200 µm have been obtained by performing a de-scan over a 7 × 7 µm 2 area with a spatial sampling δ 0 = 230 nm (961 input wave-fronts).

Local RPSF

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

I(∆ρ in , r p ) = ⟨|R in (∆ρ in , ρ out , z)| 2 W L (ρ out -ρ p )⟩ ρout (8)
where W L (ρ out -ρ p ) = 1 for |x out -x p | < L and |y out -y p | < L , and zero otherwise.

Multi-scale compensation of wave-distortions

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

L = F OV /2 q (9)
where F OV denotes the initial field-of-view.

The whole procedure is summarized in Supplementary Figure S5. At each stage of this iterative process, the starting point is the de-scanned reflection matrix R (q-1) in , obtained at the previous step, R

in being the reflection matrix recorded by our experimental set up (Fig. 1). An input distortion matrix D (q) in is deduced from R (q) in via a numerical Fourier transform (Eq. 5). A local correlation matrix of wave distortions is then built around each point r p of the field-of-view:

C in (u in , u ′ in , r p ) = D (q) in (u in , ρ out , z)D (q) * in (u ′ in , ρ out , z)W L (ρ out -ρ p ) ρout (10)
IPR is then applied to each correlation matrix C in (r p ) (see further and Supplementary Section S8). The resulting input phase laws, φin (r p ), are used to compensate for the wave distortions undergone by the incident wave-fronts:

R ′ in = F † × exp -i φin • D (q-1) in (11)
The corrected matrix R ′ in is only intermediate since phase distortions undergone by the reflected wave-fronts remain to be corrected.

To that aim, an output de-scanned matrix R ′ out (z) is deduced from the input de-scanned matrix R ′ in (z) using the following change of variable (Supplementary Figure S6):

R ′ out (ρ in , ∆ρ out , z) = R ′ in (-∆ρ out , ρ in + ∆ρ out , z) (12)
with ∆ρ out = ρ out -ρ in = -∆ρ in . An output distortion matrix is then built by applying a Fourier transform over the de-scanned coordinate:

D ′ out = R ′ out × F T ( 13
)
where the superscript T stands for matrix transpose. From D ′ out , one can build a correlation matrix C out for each point r p :

C out (u out , u ′ out , r p ) = D ′ out (ρ in , u out , z p )D ′ * out (ρ in , u ′ out , z p )W L (ρ in -ρ p ) ρ in (14)
The IPR algorithm described further is then applied to each matrix C out (r p ). The resulting output phase laws, φ out (r p ), are leveraged to compensate for the residual wave distortions undergone by the reflected wave-fronts:

R (q) out = [D ′ out exp (-iϕ in )] × F * (15)
The RPSFs displayed in Fig. 5c are extracted from the matrices R (q) out obtained at the end of each iteration of the multi-scale process. An input de-scanned matrix, combining the input and output corrections, is finally obtained by performing the following change of variables:

R (q) in (∆ρ in , ρ out , z) = R (q) out (ρ out -∆ρ in , -∆ρ in , z) (16) This matrix R (q)
in is the starting point of the next stage of the multi-scale process, and so on. Estimators of the T-matrices correspond to the cumulative function of the aberration phase laws:

T(q) in/out = T(q-1) in/out • exp iϕ (q) in/out = q k=1 exp iϕ (k) in/out (17)
Figure 5b shows the evolution of one line of the estimator T(q) out throughout the RMI process. The iterative procedure is stopped by investigating the correlation properties of this estimator (see further and Supplementary Section S13).

Iterative phase reversal algorithm.

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

φ(n) (r p ) = arg C(r p ) × exp i φ(n-1) (r p ) (18)
where φ(n) is the estimator of ϕ at the n th iteration of the phase reversal process. φ(0) is an arbitrary wave-front that initiates the process (typically a flat phase law) and φ = lim n→∞ φ(n)

is the result of IPR.

Aberration and Scattering Components of the T-matrix.

The spatial correlation of transmitted wave-fields are investigated at each depth z by computing the correlation matrix of Tout : C T = Tout × T † out . A mean correlation function Γ can be computed by performing the following average:

Γ(∆ρ, z) = ⟨C T (ρ in , ρ in + ∆ρ, z)⟩ ρ in (19)
The correlation function Γ displayed in Fig. 6a shows that the matrix Tout can be decomposed as a spatially-invariant component A out and a short-range correlated component S out . Each component can be separated by performing a singular value decomposition of Tout , such that

Tout = N p=1 s p U p V † p ( 20
)
where s p are the positive and real singular values of Tout sorted in decreasing order, U p and V p are unitary matrices whose columns correspond to the input and output singular vectors of Tout . The first eigenspace of Tout provides its spatially-invariant aberrated component:

A out = s 1 U 1 V † 1 .
The higher rank eigenstates provide the forward multiple scattering component S out . Lines or columns of the associated correlation matrix C S = Ŝout × Ŝ † out provides the isoplanatic patches displayed in Fig. 6e. this work held by CNRS (no. US11408723B2, published August 2022). All authors declare that they have no other competing interests.

45

• by a polarizer (P1) so that it is then equally reflected (sample arm) and transmitted (reference arm) by the polarized beam splitter (PBS1).

The sample beam reflected by (PBS1) is horizontally polarized. It propagates through a quarter-wave plate (QWP1), is reflected by a plane mirror (M1), whose normal axis lies along the optical axis and that is mounted on a piezoelectric actuator (PZT). The reflected beam passes again through the quarter-wave plate (QWP1). This sequence induces a polarization rotation by 90 • of the reflected beam with respect to the incident beam in the sample arm.

The reflected wave can be then transmitted through the beam splitter (PBS1) with a vertical polarization and finally focused in a secondary source plane conjugated with the source plane by means of the lens (L2) of focal length f 2 = 125 mm.

The reference beam, vertically polarized at the exit of the polarizer (P1), is transmitted by the beam splitter (PBS1), propagates through a quarter-wave plate (QWP2), is reflected by a set of two galvanometric scan mirrors, and then by the reference mirror (M2). The set of scan mirrors enables a 2D rotation of the incident wave-field by angles Θ = (Θ x , Θ y ) with respect to the optical axis. After reflection on the reference mirror (M2) and on the scan mirrors back again, the reflected beam propagates again through (QWP2). This round trip through (QWP2) enables a 90 • rotation of the polarization: Whereas the incident light is V-polarized when it enters the reference arm, it is H-polarized when exiting it (see Fig. S1a).

Therefore, the reference beam is reflected by the beam splitter (PBS1) before being focused by the lens (L2) in the secondary source plane.

Finally, in the secondary source plane, the wave-field is made of two images of the incident light orthogonally polarized and translated with respect to each other by a relative position ∆ρ = (∆x, ∆y). This lateral shift is dictated by the tilt Θ = (Θ x , Θ y ) of the reference beam: ∆x = 2f 2 tan Θ x and ∆y = 2f 2 tan Θ y . The factor 2 results from the double reflection on each scan mirror, due to the reflection on the (M2) mirror. Note also that the optical path difference between the two arms is set to zero by equalizing the length of sample and reference arms for ∆θ = 0. Regarding the horizontally-polarized beam at the exit of the lens (L3), it is reflected by the beam splitter (PBS2), passes through the quarter-wave plate (QWP3) before being focused by the microscope objective (MO1) identical to (MO2). The light is then reflected by the reference mirror (M3) placed in the focal plane of (MO2) before being collected again by the same microscope objective (MO2). The reflected light comes through the quarter-wave plate (QWP3). As in the other arm, the polarization of the reflected beam exhibits a 90 • rotation of its polarization. The beam is now vertically polarized and transmitted by the beam splitter (PBS2), before being focused on the camera with the lens (L4).

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

The CMOS camera (Adimec Quartz 2A-750, 2Mpx) records the interferogram between sample and reference beams with a spatial sampling equal to δ 0 =230 nm given the magnification M 4 . The volume of the sample from which photons can interfere with the reference beam is called the "coherence volume". Its position is dictated by the optical path difference between the reference and sample arms. Its thickness is inversely proportional to the light spectrum bandwidth 54 :

δz t = 2 ln 2 nπ λ 2 ∆λ (S1)
with λ the central wavelength of the light source and ∆λ its spectral bandwidth. In the present case, δz t ∼ 10 µm. A critical tuning of the experimental set up consists in adjusting the coherence volume with the focal plane of the microscope objective. In a volumetric sample, whose refractive index differs from that of water, the coherence volume no longer coincides with the focusing plane. This focusing defect accumulates with the transverse aberrations generated by the heterogeneities of the medium. However, it is possible to compensate for it by a fine tuning of the length of the reference arm.

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

I α (∆ρ, ρ out , z) = T 0 |e iα E out (ρ out , t) + E (ref) * out (∆ρ, ρ out , t)| 2 dt (S2)
with t the absolute time, r the position vector on the CCD screen, E out (r, τ ) the scattered wave field associated with the sample arm, E

out (r, τ ) the reference wave field; T the integration time of the CCD camera, and α an additional phase term controlled with a piezoelectric actuator placed on mirror (M1) of the first interferometer [Fig. S1a]. The interference term between the sample and reference beams is extracted from the four intensity patterns (Eq. S2) recorded at α = 0, π/2, 3π/2 and π (phase-stepping interferometry) and provides the de-scanned wave-field:

R in (∆ρ, ρ out , z) = 1 T T 0 E out (ρ out , t)E (ref) * out (∆ρ, ρ out , t)dt (S3)
Note that previous studies 25,26 reported on the passive measurement of the de-scanned reflection matrix at the surface of a scattering sample. Although the experimental set up presented in those studies shares some similarities with the current set up displayed in Fig. S1 (low-coherence interferometry), there are also several major differences. First, the two arms are illuminated by the light reflected by the sample in Ref. 25 ; hence there is no reference arm. Second, the reflection matrix is measured at the surface of the sample as a function of time-of-flight, while the current set up measures a time-gated reflection matrix as a function of depth inside the sample. At last, a much larger integration time is required to record the reflection matrix in Ref. 25 because of the absence of reference arm. In a previous work 55 , the scattering mean free path ℓ s in the cornea was measured by investigating the depth evolution of the confocal intensity. Indeed, in the single scattering regime, under the paraxial approximation and for an homogeneous reflectivity, the time-gated confocal intensity is supposed to decrease as exp(-2z/ℓ s ) if we neglect absorption losses 8,56 .

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

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

β(z) ∼ exp [-2z/(3ℓ s )] (S4)
In the present case, this confocal ratio β has been measured as follows:

β(z) = 1 - min ∆ρ {I(∆ρ, z)} max ∆ρ {I(∆ρ, z)} (S5)
This estimator β relies on the fact that the multiple scattering component of the RPSF exhibits a flat background such that it can estimated with the minimum of I(∆ρ, z). This hypothesis is wrong at shallow depth since the diffuse halo grows as √ Dt. Nevertheless, beyond ℓ s or so (here 100 µm), the multiple scattering background can be considered as flat as illustrated by Fig. 3d of the accompanying paper.

S3. THEORETICAL EXPRESSION OF THE DE-SCANNED MATRIX

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

The wave field E out (ρ out , z) reflected by the sample arm in the camera plane can then be expressed as follows 29 : out (ρ out , z) 29 :

E out (ρ out , z, ω) = Σ 0 Σρ H(ρ out , ρ s , z)γ(ρ s , z)H(ρ s , ρ 0 , z)E 0 (ρ 0 , ω)dρ s dρ 0 . ( S6
E (ref) out (ρ out , ∆ρ, z) = Σ 0 H ref (ρ out -ρ 0 )E 0 (ρ 0 + ∆ρ)dρ 0 . (S7)
where H ref is the impulse response associated with the reference arm (way and return path) that we assume as spatially

-invariant [H ref (ρ out , ρ 0 ) = H ref (ρ out -ρ 0 )].
The de-scanned wave-field is obtained by extracting the interference term between the reflected wave-fields coming from the sample and reference arms:

R in (ρ out , ∆ρ, z) = ⟨E out (ρ out , ω)E (ref) * out (ρ out , ω)⟩ (S8)
Assuming a spatially-incoherent incident wave-field [⟨E 0 (ρ 0 )E * 0 (ρ ′ 0 )⟩ = I 0 δ(ρ 0 -ρ ′ 0 )] and injecting Eqs. S6 and S7 into the last equation leads to the following expression for the coefficients of R in :

R in (ρ out , ∆ρ, z) = I 0 H out (ρ out , ρ, z)γ(ρ s , z)H in (ρ s , ρ out + ∆ρ, z)dρ (S9)
with For ∆ρ = 0 (conventional FFOCT set up), the recorded wave-field (Eq. S9) is equivalent to a time-gated confocal image 29 . It can actually be expressed as the convolution between the sample reflectivity γ and the confocal PSF

H out ≡ H and H in ≡ H ⊙ H * ref . ( S10
H in × H out : R in (ρ out , ∆ρ = 0, z) = I 0 H out (ρ out , ρ s , z)H in (ρ s , ρ out , z)γ(ρ s , z)dρ s . (S11)
On the one hand, the confocal nature of the recorded wave-field implies a transverse resolution δρ 0 ∼ λ/4NA. On the other hand, the axial resolution is either controlled by the thickness δz t of the coherence volume or the depth-of-field δz 0 of the microscope objective: δz 0 = nλ/NA 2 .

In the present case, δz 0 ∼ 1 µm< δz t ∼ 10 µm. The axial resolution is thus given by the depth-of-field. δρ 0 and δz 0 thus dictate the values of the transverse and axial sampling of the de-scanned matrix R in in our experiment.

S4. RELATION BETWEEN THE DE-SCANNED MATRIX AND THE FOCUSED REFLECTION MATRIX

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

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

R(ρ out , ρ in , z) = I 0 H(ρ out , ρ s , z)γ(ρ s , z)H(ρ s , ρ in , z)dρ (S12)
A strict equality between Eqs. S9 and S12 is only obtained if H in ≡ H. This condition is fulfilled only for a perfect reference arm:

T ref ≡ 1 k<N A and H ⊙ H * ref ≡ H.
In theory, the incoherent set up of Fig. S1a is thus equivalent to the fictitious coherent set up of Fig. S1b.

R in (∆ρ, ρ out ) = R(ρ out + ∆ρ, ρ out , z) (S13)
In reality, the reference arm always exhibits aberrations such as a slight defocus of the reference mirror M3 in Fig. S1b or a slight defocus of the reference beam in the secondary source plane at the output of first interferometer.

S5. THE DISTORTION MATRIX

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

D in (z) = F × R in (z), (S14)
or in terms of matrix coefficients,

D(u in , ρ out , z) = ∆ρ R in (∆ρ, ρ out , z) exp -i 2π λf u in • ∆ρ . (S15)
Injecting Eq. S9 into the last equation yields

D(u in , ρ out , z) = ∆ρ I 0 H out (ρ out , ρ s , z)γ(ρ s , z)H in (ρ s , ρ out +∆ρ, z) exp -i 2π λf u in • ∆ρ dρ s . (S16)
In a previous paper 19 , we showed that a singular value decomposition of D enables to decompose the field-of-view into isoplanatic modes and extract the associated aberration phase laws. However, this demonstration was based on the condition that the D-matrix is dominated by its correlations in the focal plane. This is the case for a specular reflector such as a resolution target or a medium of continuous reflectivity but no longer valid for a random distribution of heterogeneities like in the opaque cornea under study. In the accompanying paper, we propose a more general solution to overcome aberrations and scattering in optical microscopy: An iterative multi-scale analysis of wave distortions.

To that aim, the field-of-detection should be subdivided into overlapping regions that are defined by their central midpoint r p = (ρ p , z p ) and their spatial extension L. All of the distorted components associated with focusing points ρ out located within each region are extracted and stored in a local distortion matrix D ′ in (r p ):

D ′ (u in , ρ out , r p ) = D(u in , ρ out , z p ) W L (ρ out -ρ p ), ( S17
)
where W L (x, y) = 1 for |x| < L and |y| < L, and zero otherwise.

At this stage, a local isoplanatic assumption shall be made over each region of size L.

This hypothesis implies that the PSFs H in/out are invariant by translation in each region.

This leads us to define local spatially-invariant PSFs

H (l)
in/out around each central midpoint r p such that:

H in/out (ρ s , ρ in/out , z p ) = H (l) in/out (ρ s -ρ in/out , r p ). ( S18
)
Under this assumption, Eq.S16 can be rewritten as follows:

D ′ (u in , ρ out , r p ) = T in (u in , r p ) transmittance γ(ρ s + ρ out , z)H (l) out (ρ s , r p ) exp -i 2π λf u in • ρ dρ s virtual source . (S19)
Around each point r p , the aberrations can be modelled by a transmittance T in (u in , r p ). This transmittance is the Fourier transform of the input PSF H (l) in (ρ s , r p ):

T in (u in , r p ) = H (l) in (ρ s , r p ) exp 2π λf u in • ρ s dρ s (S20)
The physical meaning of this last equation is the following: Each distorted wave-field corresponds to the diffraction of a virtual source synthesized inside the medium modulated by the transmittance T in of the sample between the focal and pupil planes. Each virtual source is spatially incoherent due to the random reflectivity of the medium, and its size is governed by the spatial extension of the output focal spot. The idea is now to smartly combine each virtual source to generate a coherent guide star and estimate T in independently from the sample reflectivity.

S6. COVARIANCE MATRIX OF WAVE DISTORTIONS

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

C in (u in , u ′ in , r p ) = N -1 W ρout D ′ (u in , ρ out , r p )D ′ * (u ′ in , ρ out , r p ) (S21)
The matrix C in (r p ) can be decomposed as the sum of its ensemble average, the covariance matrix ⟨C in ⟩ (r p ), and a perturbation term δC in (r p ):

C in (r p ) = ⟨C in ⟩ (r p ) + δC in (r p ). ( S22
)
The intensity of the perturbation term scales as the inverse of the number N W = (L/δρ 0 ) 2 of resolution cells in each sub-region 38,58 :

|δC in (u, u ′ , r p )| 2 = |C in (u, u, r p )| 2 N W (S23)
This perturbation term can thus be reduced by increasing the size L of the spatial window W L , but at the cost of a resolution loss.

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

⟨γ(ρ s , z)γ * (ρ ′ s , z)⟩ = ⟨|γ| 2 ⟩δ(ρ s -ρ ′ s ), (S24)
with δ, the Dirac distribution, the covariance matrix can be expressed as follows 21 :

⟨C in ⟩ (r p ) = [T in (r p ) • F] × C H (r p ) × [F • T in (r p )] † , (S25)
or in terms of matrix coefficients,

⟨C in ⟩ (u in , u ′ in , r p ) = T in (u in , r p )T * in (u ′ in , r p ) dρ s H (l) out (ρ s , r p ) 2 exp -i 2π λf (u in -u ′ in ) • ρ s = T in (u in , r p )T * in (u ′ in , r p ) T out u ⊛ T out (u in -u ′ in , r p ) =C H (u in ,u ′ in ,rp) . ( S26
)
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 (l) out (ρ s , r p )| 2 .
The covariance matrix ⟨C in ⟩ (r p ) thus corresponds to the same experimental situation but for a virtual reflector embedded into the heterogeneous cornea under study.

S7. ITERATIVE TIME REVERSAL

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

ψ + u = C in × ψ - u (S27)
and

ψ - u = ϵψ + * u . (S28)
with ϵ the reflectivity of the phase conjugating mirror. Combining the two previous equations leads to the following eigenequation:

ψ - u = |ϵ| 2 C * in C in × ψ - u . (S29)
The ITR process has thus eigenmodes which can be determined by the diagonalization of the time reversal operator C * in C in . In particular, the first eigenvector U in of C * in C in , which is also the first singular vector of D in , corresponds to the wave-front that optimizes the energy backscattered by the virtual reflector.

If the virtual reflector was point-like, this wave-front would be a perfect estimator of T in .

Its phase conjugate would perfectly compensate for aberrations and focuses through the heterogeneous medium onto the point-like target 59,60 . However, here the virtual guide star is enlarged compared to the diffraction limit. This wave-front is of finite angular support δu c and tends to focus on the virtual reflector but with a resolution width δρ c ∼ λf /δu c larger than the diffraction limit 38 (see Fig. S3). Its phase is thus a good estimator of T in over the angular domain δu c but absolutely not elsewhere. This assertion is illustrated by Figs. S3b andc that show the modulus and phase of the first singular vector U in of D in at depth z = 200 µm in the cornea. As anticipated, the modulus of U in exhibits a main central lobe at small spatial frequencies (delimited by a circle white line in Fig. S3b) but is extremely low at high angles of incidence. This means that the phase of U in is a good estimator for |u in | < δu c but is not reliable beyond u c (Fig. S3c). matrix C in but a similar demonstration can be performed at output. For sake of lighter notation, the dependence over r p is omitted in the following.

To understand the parameters controlling the error δT in between Tin and T in , one can express Tin as follows:

Tin = exp jarg C in × Tin = C in × Tin ||C in × Tin || (S32)
By injecting Eq. S22 into the last expression, Tin can be expressed, at first order, as the sum of its expected value T in and a perturbation term δ Tin :

Tin = ⟨C in ⟩ × T in ||⟨C in ⟩ × T in || =T in + δC in × T in ||⟨C in ⟩ × T in || ≃δ Tin . ( S33
)
The bias intensity can be expressed as follows:

|δ Tin | 2 = T † in × δC † in × δC in × T in T † in × ⟨C in ⟩ † × ⟨C in ⟩ × T in (S34)
Using Eq. S23, the numerator of the previous equation can be expressed as follows:

T † in × δC † in × δC in × T in = M 2 ⟨|δC in (u, u ′ )| 2 ⟩ = M 2 ⟨|⟨C in ⟩(u, u)| 2 ⟩/N W . (S35)
Injecting Eq. S26 into the last equation leads to the following expression for the numerator of Eq. S34:

T † in × δC † in × δC in × T in = M 2 T out u ⊛ T out (0) 2 /N W . (S36)
The denominator of Eq. S34 can be expressed as follows:

T † in × ⟨C in ⟩ † × ⟨C in ⟩ × T in = M 2 u T out uout ⊛ T out (u) 2 (S37)
The bias intensity is thus given by:

|δT in | 2 = T out uout ⊛ T out (0) 2 N W u T out u ⊛ T out (u) 2 (S38)
In the last expression, we recognize the ratio between the coherent intensity (energy deposited exactly at focus) and the mean incoherent intensity. This quantity is known as the coherence factor in ultrasound imaging 39,58 :

C out = u T out u ⊛ T out (u) T out u ⊛ T out (0) = |H (l) out (ρ = 0)| 2 δ -2 R dρ|H (l) out (ρ)| 2 (S39)
In the speckle regime (Eq. S24) and for 3D imaging, the coherence factor C ranges from 0, for strong aberrations and/or multiple scattering background, to 4/9 in the ideal case 40 .The bias intensity can thus be rewritten as:

|δT in (u)| 2 = 1 C 2 out N W (S40)
This last expression justifies the multi-scale analysis proposed in the accompanying paper.

A gradual increase of the focusing quality, quantified by C, is required to address smaller spatial windows that scale as N W . Following this scheme, the bias made of our T-matrix estimator can be minimized and the iterative phase reversal algorithm converges towards a satisfying estimator.

S10. NUMERICAL VALIDATION OF THE ITERATIVE PHASE REVERSAL

PROCESS

The IPR algorithm is now validated by means of a numerical simulation. The numerical simulation emulates an imaging experiment in an epi-detection configuration, as depicted in matrix R, as described in section S5, and the estimations of the input and output aberration transmittances Tin/out (u in/out ) are computed using the IPR process described in the Methods section of the accompanying paper.

Figures S4c andd show the estimated transmittances Tin/out , respectively, when the whole FOV is considered. A strong similarity is observed with the ground truth up to a phase ramp 38 (Fig. S4b). The corresponding RSPFs after each correction are displayed in Figs. S4f andg. A diffraction-limited resolution is obtained at the end of the RMI process, which validates the IPR algorithm.

One can go further by investigating the convergence of the process as a function of N W , the size of the spatial window considered for the computation of the correlation matrices C in/out . The similarity between the estimators Tin/out and the ground truth T is evaluated by the normalized scalar product Q in/out = N -1 u Tin/out × T † , or, in terms of matrix coefficients.

Q in/out = N -1 u u Tin/out (u)T ⋆ (u). ( S43
)
The evolution of Q in and P out is displayed as a function of N W in Fig. S4h. The convergence can be considered as fulfilled for Q > 0.9, i.e N W ∼ 500, which is roughly the number of resolution cells contained in the final spatial windows (L = 6 µm) in our experiment.

This convergence rate is directly related to the bias of Tin/out (Eq. S40). To show it, let us first express the intensity bias |δT (u)| 2 as a function of the phase error δϕ(u) exhibited by the estimator T (u) with respect to T (u), following the same formalism as in section S9:

|δT (u)| 2 = |1 -exp [iδϕ(u)] | 2 ∼ δϕ≪1 |1 -(1 + δϕ)| 2 ∼ δϕ≪1 |δϕ| 2 (S44)
On the other hand, the scalar product as a function of the phase error δϕ writes as such:

Q in/out = N -1 u u exp [iδϕ(u)] (S45)
The sum over the points in the Fourier plane u can be replaced by an ensemble average, since N >> 1:

Q in/out = ⟨exp [iδϕ(u)]⟩ (S46)
Assuming a small phase error (δϕ ≪ 1),

Q in/out ∼ 1 + i ⟨δϕ⟩ - ⟨(δϕ) 2 ⟩ 2 (S47) ∼ 1 - ⟨(δϕ) 2 ⟩ 2 (S48)
since ⟨δϕ⟩ = 0. Combining this last expression with Eq. S44 leads to:

Q in/out ∼ 1 - ⟨|δT | 2 ⟩ 2 (S49)
According to Eq. S9, 1 -Q in/out should therefore scale as the inverse of the number of independent resolution cells contained in the spatial window: 1 -P in/out ∝ N -1 W . To highlight this scaling law, 1 -Q in/out can be plotted in log-log scale as a function of N W (Fig. S4h).

A slope p close to 1 is obtained both at input and output: p in = -0.95 and p out = -1.1, confirming that the bias on the aberration estimation scales with the inverse of the number of independent resolution cells in the field of view (Eq. S9). Another interesting observation is the lower bias observed at output in Fig. S4h,i. Indeed, the first correction at input increases the coherence factor C in and reduces the size of the virtual guide star when investigating wave distortions at output. This gain in focusing quality improves the sharpness of the estimator T, as already highlighted by the scaling of |δT | 2 as the inverse square of the coherence factor in Eq. S9.

In the present numerical simulation, the isoplanicity assumption makes the IPR algorithm converging towards an appropriate solution in one iteration at input and output. In the experiment, the situation is more complex since aberrations are spatially-distributed. In that case, an iterative compensation of wave distortions aver a multiple scale is required. The corresponding strategy is explained in the next Section. Eq. 17

Eq. 17 The multi-scale compensation of wave distortions consists in dividing by two the lateral extension L of the spatial windows W at each step. The full process is described in the Methods section of the accompanying paper and summarized in a flowchart displayed in Fig. S5. At each step, the correction process is iterated both at input and output of the reflection matrix (left and right parts of Fig. S5). Mathematically, the transfer between the input and output de-scanned bases is performed by a change of variable (Eqs. 12 and 16) illustrated by Fig. S6. In particular, Fig. S6f shows that the output de-scan matrix R out cannot be fully retrieved. A set of coefficients cannot be determined in its corners and are arbitrarily fixed to zero. They correspond to de-scanned coordinates (ρ in , ∆ρ ′ ) associated with points ρ out = ρ in + ∆ρ ′ outside of the initial field-of-detection. To avoid the potential detrimental impact of such zero coefficients on the estimation of the T-matrix, the output correlation matrix C out is only computed over points ρ in that are associated with a full de-scan wave-field, i.e points ρ in such that R out (ρ in , ∆ρ ′ ) ̸ = 0 for each de-scan position ∆ρ ′ . In a previous work 19 , the compensation of wave-front distortions was performed in one single step and on a single side (output). The low spatial sampling of the reflection matrix at input explained this minimalist strategy. In the accompanying paper, the de-scanned measurement of the reflection matrix provides the same sampling of the wave-field at input and output. An alternate compensation of wave distortions is therefore possible and actually critical if one wants to converge towards a sharp estimator of the T-matrix. Indeed, as shown by Eq. S40, the bias of this estimator on one side (input/output) directly depends on the focusing quality on the other side (output/input) since it controls the blurring of the virtual andf) and that RMI tends to realign in phase by means of a digital optical phase conjugation.

Figure S9b clearly shows that the confocal gain increases with depth z. Indeed, multiple scattering becomes predominant in depth and the transmission phase laws become more and more complex. Note, however, that given the complexity of phase laws displayed in Figs. 4e andf, we could have expected a larger confocal intensity enhancement. This moderate gain reference arm is isoplanatic, it can be extracted by considering the first eigenstate of Tin (see Supplementary Section S15). The phase of the pupil singular vector U 1 displayed in Fig. S10 is an estimator of T ref . Not surprisingly, it mainly corresponds to a spherical aberration phase law. One can subtract this reference phase to Tin in order to build a matrix T′ in = Tin • U * 1 . One can expect the scalar product P ′ (ρ p , z) between T′ in and Tout to be increased compared to its initial value P (ρ p , z) (see comparison between Figs. S10a andc). This is actually what we observe even though the scalar product P ′ remains smaller than 0.7 (Fig. S10c). It means that the spherical aberration law induced by the reference arm account partially for the mismatch between Tin and Tout .

The residual mismatch can be explained by the fact that the aberration induced by the reference arm is not strictly isoplanatic. Misalignment between sample and reference arms manifests as a transverse shift of the RPSF that varies across the field-of-view as illustrated by Fig. 4c. Field curvature can also induce spatially-varying aberrations that our approach can address but they are difficult to discriminate from sample arm aberrations. Last but not least, another phenomenon that can contribute to this discrepancy between input and output aberration phase laws is the bias of our T-matrix estimator, especially for small spatial windows L as explained in Supplementary Section S9.

S15. ABERRATION AND SCATTERING COMPONENTS OF THE T-MATRIX

In a previous work, Badon et al. showed how the singular value decomposition (SVD) of the D-matrix provided a decomposition of the field-of-view into isoplanatic modes in the case of a specular object. In the present paper, this property does not hold since we cope with a random distribution of heterogeneities. In this regime, this is the SVD of the T-matrix estimator that enables a mapping of isoplanatic modes. As recently shown in an ultrasound study 38 , the complexity of the associated aberration phase laws increases with the rank of the corresponding singular values while the spatial extension of the isoplanatic mode decreases. As we will show below, this complexity can be quantified by a vorticity degree of the associated transmittance, quantity that has a direct link with the occurrence of multiple scattering paths involved in the trajectory of the wave from the focal plane to the camera sensors. , such that ∇ • ∇ϕ (rot) k = 0, induced by forward multiple scattering trajectories. Indeed, this curl component is a manifestation of optical vortices that necessarily originate from, at least, three interfering beams and thus suppose several optical trajectories, hence multiple scattering. A degree α of vorticity can be assessed by looking at the ratio between the energy of each contribution, such that:

α k = |∇ ∧ ∇ϕ k | 2 |∆ϕ k | 2 (S57)
This degree of vorticity is displayed for each eigenstate in Fig. S11b. Although it shows some fluctuations, α k tends to increase with the rank k of eigenstate. It thus seems to indicate that the higher order eigenstates associated with smaller isoplanatic patches also exhibit a higher degree of vorticity, which is a manifestation of forward multiple scattering paths.

S16. COMPARISON BETWEEN MATRIX IMAGING AND ADAPTIVE OPTICS

The discrimination between aberration and scattering phenomena can be also applied directly to the estimated T-matrix. The interest of such a decomposition is to show the superiority of matrix imaging with respect to conventional AO. Indeed, the latter approach generally relies on Shack Hartmann sensors that only give access to the phase gradient of reflected wave-fronts. Standard numerical integration of this quantity gives access to the irrotational (i.e aberrated) component of the wave-front but generally not to the scattering components of the wave-front that exhibits a wealth of optical vortices 42 . On the contrary, the interferometric measurement of the reflected wave-field gives access to this scattering component which is crucial for deep imaging. components. The access to the latter component is decisive for the compensation of wave distortions since it greatly contributes to the improvement of the confocal image (Fig. S12b).

This can be quantified by the confocal gain exhibited at the end of the matrix imaging process (Fig. S12c) and the corresponding RPSF (Fig. S12d). While the access to the curl component of the focusing laws allows us to reach a confocal gain up to 13 dB Fig. S12c 1 ), conventional AO would only allow a compensation of low-order aberrations, giving rise to a weak confocal gain (< 3 dB, Fig. S12c 2 ). effects (Fig. S13b 3 ), which is a manifestation of a lack of correlations of the estimator T between adjacent windows. On the contrary, the IPR leads to an estimator whose scattered component Ŝ exhibits the expected memory effect for this size of spatial windows (Fig. S8d).

Figure S13 thus shows the superiority of IPR in terms of spatial resolution for the T-matrix estimation.

S18. COMPARISON BETWEEN CLASS AND RMI

For the last five years, an alternative matrix approach has been developed by W. Choi and colleagues, this is the so-called CLASS algorithm 20,23,24 . Based on the recording of a time-gated reflection matrix, it exploits its input-output correlations in the plane wave basis image displays a more homogeneous reflectivity across the field-of-view. This difference can be understood by the maximization of the confocal energy operated by CLASS which gives more weight to the most echogenic scatterers to the detriment of weaker reflectivity regions. Another difference between CLASS and IPR approach lies in the higher resolution capability of the latter method. This assertion is supported by Fig. S15 that compares the results of CLASS and IPR for a smaller spatial window (L = 5.75 µm). A local RPSF obtained via CLASS is displayed in Fig. S15b 2 . It shows an important background, a manifestation of an imperfect compensation of forward multiple scattering. On the contrary, the RPSF obtained at the same location via IPR is close to be ideal (Fig. S15c 2 ). When looking into details to the CLASS image (Fig. S15b 1 ), strong vignetting effects can be observed while the IPR image exhibits a continuous reflectivity (Fig. S15c 1 ).

This result can be understood by comparing the spatial correlation properties of the scattering component S of the T-matrix obtained by CLASS and IPR (see Methods of the accompanying paper). Figure S16 shows this correlation map for the mid-point r p of the area displayed with a white square in Fig. S15. While the S-matrix derived by IPR preserves a short-range correlation between neighbouring windows (see Fig. S16c), the CLASS algorithm leads to a fully spatially incoherent estimator Ŝ (see Fig. S16a). This observable clearly shows that the IPR estimator leads to a coherent (i.e physical) compensation of multiple scattering while CLASS leads to an incoherent correction (i.e bucket-like). The result of

FIG. 1 .FIG. 2 .

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

FIG. 3 .

 3 FIG. 3. Different Stages of Matrix Imaging. a. The focused R-matrix contains the set of impulse responses R(ρ in , ρ out , z) between an array of point sources ρ in and detectors ρ out lying in planes conjugated with the focal plane of the microscope objective. b. The interferometric set up displayed in Fig. 1 allows a de-scanned measurement of R by scanning the relative position ∆ρ in = ρ in -ρ out . c. Each column of the recorded matrix R in (Eq. 1) corresponds to the RPSF measured by each camera pixel. d. A spatial Fourier transform over ∆ρ in provides the distortion matrix D in (Eq. 5) linking each camera pixel with wave-front distortions seen from the input pupil plane (u in ). e. The correlation matrix C in between those wave-fronts mimics the time-reversal

FIG. 4 .

 4 FIG. 4. Time-gated Transmission matrix for Local Compensation of Forward Multiple Scattering. a,b. Confocal field of view before and after the correction process at 200 µm-depth, respectively (scale bar: 50 µm). c,d. Maps of the local reflection point-spread functions (RPSFs) (de-scan field-of-view: 7 × 7 µm 2 ) over the field of view, before and after the correction process, respectively. e,f. Sub-part of matrices, Tin and Tout , respectively, for the area delimited by the square box in panels a-d.

FIG. 5 .

 5 FIG. 5. Multi-scale analysis of wave distortions. a. The entire field-of-view is 138 × 138 µm 2 . At each step, it is divided into a set of spatial windows whose dimension gradually decreases: from 138, 100, 50, 25, 13 to 6 µm b. Evolution of the transmitted wave-front T (u out , r p ) for one point r p of the field-of-view at each iteration step. c. Corresponding local RPSF at r p before and after compensation of aberration and scattering using digital phase-conjugation of the optical transfer function displayed in panel b (scale bar: 2 µm). Data are from the cross-section at 200 µm depth within the sample.

FIG. 6 .

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

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

  After the Michelson interferometer, the two orthogonally polarized twin beams enter a Michelson interferometer with two identical microscope objectives in both arms (a configuration known as a Linnik interferometer) [Fig.S1b]. They are again collimated by a lens (L3) of focal length f 3 =200 mm. The two lenses (L2) and (L3) thus constitute a 4f system which compensates the effects of diffraction between the two interferometers.The vertically polarized light (sample beam) is transmitted by a polarized beam splitter cube (PBS2), propagates through a quarter-wave plate (QWP4) before being focused in the focal plane of an immersion microscope objective (MO2, Nikon, 60×, NA=1.0). The light reflected by the sample is then collected by (MO2) and propagates again through the quarter-wave plate (QWP4). Because single scattering tends to preserve polarization, the corresponding wave-field undergoes a 90 • polarization rotation and gets reflected by the beam splitter (PBS2) before being focused in the plane of the camera using the converging lens (L4) of focal length f 4 =200 mm. The combination of this lens (L4) with the microscope objective (MO1) entails a magnification M 4 of 60.

  FIG. S1. De-scanned measurement of the reflection matrix. a. Experimental set up. P: polarizer, L: lens, QWP: quarter-wave plate, M: mirror, PZT: piezo-electric actuator, PBS: polarisation beam splitter, MO: microscope objective. The apparatus is made up of two parts (a 1 ,a 2 ). a 1 . Michelson interferometer illuminated by an incoherent light source at its input and generating two twin incoherent beams of orthogonal polarization and laterally shifted from each other at its output. The polarised beam splitter (PBS1) separates the impinging light into a reference path (in blue) and a sample path (in red). The scan mirrors (a 3 ) tilt the reference beam by angles ∆Θ x and ∆Θ y in both transverse directions x and y. This allows the scan of the point-spread function along the de-scanned coordinates ∆ x and ∆ y in a plane conjugate to the reference mirror plane. a 2 . Michelson interferometer with microscope objectives (MO) in both arms (Linnik configuration). Both beams have orthogonal polarizations and each interferometer arm includes a quarter-wave plate (QWP). The output beams are collected by the L4 lens and interfere on the camera after having been projected on a 45 • -rotated polarizer (P2). b. Equivalent layout in the case of a coherent measurement. The source plane, the focal plane, and the camera planes are conjugated. Displacing a point source ρ in in the source plane discretely scans the focal plane inside the sample. The illuminated area is imaged in the camera plane; in an epi-detection configuration.

Figure

  Figure S2 displays the depth evolution of the estimator β(z). It exhibits an exponential decay in the stroma beyond z = 100 µm. The decay rate decreases beyond z = 170 µm because our estimator of β(z) starts to be impacted by the experimental noise [see Fig. 3d of the accompanying paper]. Therefore, the fit of β(z) with Eq. S4 is performed from z = 100 to z = 170 µm. We find ℓ s ∼ 35 µm.

) E 0

 0 (ρ 0 , ω) is the incident wave-field in the secondary source plane Σ 0 at frequency ω. Light propagation between Σ 0 and the focal plane Σ ρ is described by the impulse response H(ρ 0 , ρ s , z) between a point in the secondary source plane at transverse coordinate ρ 0 and a point at transverse coordinate ρ s in the focal plane and at depth z inside the sample. It accounts for sample-induced aberrations. γ(ρ s , z) represents the sample reflectivity at depth z. By spatial reciprocity, the propagation of the reflected wave-field from the sample to the detector plane is also modelled by the impulse response H(ρ s , ρ out ). The relatively narrow bandwidth (∆λ ≪ λ) of the light source and the use of achromatic optical elements (lens, beam splitter, quarter wave plate) allow us to neglect the dependence of H on frequency ω.Replacing γ(ρ s , z) by a uniform reflectivity in Eq. S6 and taking into account the lateral shift of the reference wave-field induced by the galvanometer M2 [Fig.S1] leads to the following previous expression for E(ref)

)

  The symbol ⊙ stands for the convolution product over the variable ρ in . The last equation means that: (i) the output focusing matrix, H out , and the associated T-matrix, T out = F × H out , only grasp the sample-induced aberrations; (ii) the input focusing matrix, H in , and the associated T-matrix, T in = F × H in , also contain the aberrations undergone by the incident and reflected reference beams (Supplementary Section S14).

  FIG. S3. Iterative Time Reversal vs. Iterative Phase Reversal. (a) The first eigenstate of C in corresponds to the eigenmode that would arise between a phase conjugating mirror in the pupil plane and the virtual reflector. (b,c) Absolute value and phase of U 1 the first eigenvector of C in at the first iteration of the multi-scale analysis (z = 200 µm). (d) The IPR process converges towards the wave-front Tin that would be obtained if an iterative phase reversal mirror was used to focus on the virtual reflector.(e,f ) Absolute value of C in × Tin and phase of Tin at the first iteration of the multi-scale analysis (z = 200 µm).

Fig. S4a .

 S4a Fig.S4a. The experimental conditions (numerical aperture, focal length, 'etc.) are identical to our experiment. The field-of-view contains N = 61 × 61 independent resolution cells. For sake of simplicity, a plane object of random complex reflectivity γ(ρ) is considered in the focal plane of the microscope objective and the isoplanatic assumption is also made. Under these assumptions, the coefficients of the reflection matrix can be expressed in the pupil basis as follows:

FIG. S5 .

 S5 FIG. S5. Flowchart of the multi-scale matrix imaging process.

  FIG. S6. Conversion of the reflection matrix from the input to the output de-scanned basis. a Reflection matrix R in in the input de-scanned basis. b Each column reshaped in 2D corresponds to the input focal spot de-scanned with respect to the output focusing point ρ out . c Each focal spot can be re-expressed in the laboratory frame (ρ in ) and d stored in the canonical reflection matrix R. e Each line of R can be reshaped in 2D corresponds to the output focal spot in the focused basis (ρ out ). f Each output focal spot can be de-scanned with respect to the input focusing point ρ in and stored in the output de-scan matrix R out .

S13.

  FIG. S9. Confocal gain provided by the matrix imaging process. ab. Transverse crosssection of the confocal gain observed for the en-face images displayed in Fig. 3b at depths 50 µm and 250 µm within the cornea [scale bar: 50 µm]. c. Longitudinal cross-section of the confocal gain observed by comparing the B-scan displayed in Fig.3f with its original version shown in Fig. 3c. In each panel, the color scale is in dB.

Figure

  Figure S9 shows the enhancement of the confocal peak before and after RMI. It reaches a maximal value of 30. This gain should scale, in amplitude, as the number P c of independent coherence grains exhibited by the T-matrix in the pupil plane (see, for instance, Figs. 4e

4 FIG

 4 FIG. S11. Singular value decomposition of the T-matrix at depth z = 50 µm. a Normalized singular values of Tout . b Degree of vorticity associated with each pupil singular vector U k of Tout . c Four first eigenstates of Tout : Phase of the transmittance U k and modulus of the corresponding isoplanatic modes V k .

Figure

  Figure S12 illustrates this assertion by first showing the decomposition of the input and output phase laws (Figs. S12a 1 ) into their irrotational (Figs. S12a 2 ) and curl (Figs. S12a 3 )

ac

  FIG. S12. Decomposition of the focusing law ϕ extracted by RMI at depth z = 200 µm.. a Examples of input and output phase laws (1) extracted by RMI and their decomposition into irrotational (2) and curl (a3) components. These phase laws correspond to the center of the white rectangle in b 1 . b Confocal image obtained using the corresponding focusing laws [scale 50 µm]. c Associated confocal gain. c Example of local RPSF obtained in the white rectangle displayed in panel b 1 ).

  FIG. S13. Multi-scale compensation of wave distortions: Iterative time reversal vs. iterative phase reversal. a Ratio between confocal gains obtained at z = 200 µm using a ITR or IPR process. Left colorbar: Ratio between IPR and SVD confocal gains. Right colorbar: Ratio between SVD and IPR confocal gains. b,c Confocal image obtained via IPR and ITR, respectively. The top (1) and bottom (2) figures correspond to size of spatial windows L = 25 and 6 µm, respectively.

FIG. S15 .

 S15 FIG. S15. Comparison between CLASS and RMI for small isoplanatic patches. a. Original confocal image (a 1 ) and mean RPSF (a 2 ) at z =200 µm. b, c. Corresponding CLASS (b 1 ) and IPR (c 1 ) images with their mean RPSFs (b 2 and c 2 , respectively). The size of spatial windows is L = 6 µm.

Fig. S16 accounts

  Fig. S16 accounts for the vignetting effects observed for the CLASS image in Fig. S15b 1 .

  H in (ρ s , ρ in , z) and H out (ρ s , ρ out , z) are the local input and output point spread functions (PSFs) at points (ρ in , z) and (ρ out , z) , respectively. This last equation confirms that the central line of R in (∆ρ

in = 0), i.e. the FFOCT image, results from a convolution between the sample reflectivity γ and the local confocal PSF, H in × H out .

S11. MULTI-SCALE COMPENSATION OF WAVE DISTORTIONS

			change of variables
		R in FT Eq. 5 (q) D in spatial window (q)	e ect no memory Eq. 16 STOP q+1 ← q iteration	T out (q) R out	Eq. 15 conjugation digital phase
		Eq. 10		φ out (r p ) (q)
		C in (r p )	iterative phase reversal
				Eq. 18
	iterative phase reversal
		Eq. 18		C out (r p )
		(q) φ in (r p )		spatial window Eq. 14
	digital
	phase conjugation	T in			out D'
	Eq. 11			FT Eq. 13
		R' in	change of variables Eq. 12	R' out

TABLE S1 .

 S1 R in (∆ρ in/out , ρ out/in , z)] de-scanned matrix at input/outputD in/out = [D(u in/out , ρ out/in ), z]input/ouput distortion matrixH in/out = [H in (ρ s , ρ in/out , z)] input/output focusing matrix F = [F (u, ∆ρ)]Fourier transform operatorC in/out = [C in/out (u in/out , u ′ in/out , r p )] input/output pupil correlation matrix of D in C H = [C H (u, u ′ )] correlation matrixassociated with a virtual reflector of reflectivity |H(ρ s )| 2 T in/out = [T in/out (u in/out , r p )] input/output transmission matrix Tin/out = [T in/out (u in/out , r p )] estimator of the input/output transmission matrix ϕ in/out = [ϕ in/out (u in/out , r p )] phase of the input/output transmission matrix estimator A in/out = [T in/out (u in/out , r p )] input/output aberration matrix S in/out = [S in/out (u in/out , r p )] input/output forward multiple scattering matrix Glossary

of the operators used in this study

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

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

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

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

Supplementary Information

This document provides further information on: (i) the experimental set up; (ii) the measurement of the scattering mean free path; (iii) the theoretical expression of the descanned matrix; (iv) The light source is a broadband LED (Thorlabs M850LP1, λ • = 850 nm, ∆λ = 35 nm), is placed in a plan conjugated with the focal plane within the sample, so as to illuminate the focal plane with the image of the source. The source's illumination pattern is not uniform and is smaller than the maximal extension of the field of view allowed by the microscope objectives' numerical aperture. As a result, the sample's intensity in the focal plane is modulated by the image of the source. In order to get a uniform illumination of the whole field of view, we set up a pseudo-Koehler illumination apparatus: An aspheric lens and a diaphragm are placed right in front of the source, such that the incident beam is collimated in the diaphragm plane. This plane is considered the source plane, and is conjugated to the sample plane by (L1). This way, the image of the source is defocused in the sample plane.

This ensures an incoherent [START_REF] Mertz | Introduction to Optical Microscopy[END_REF] , yet uniform, illumination of the field of view.

The incident light is collimated using a converging lens (L1) with a focal length f 1 = 150 mm. The beam transmitted through this lens (L1) is linearly polarized at

S8. ITERATIVE PHASE REVERSAL

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

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

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

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

S9. BIAS OF THE T-MATRIX ESTIMATOR

The IPR process assumes the convergence of the correlation matrix C in (Eq. 10) towards its ensemble average ⟨C in ⟩, the covariance matrix 21,38 . In fact, this convergence is never fully realized and C in should be decomposed as the sum of this covariance matrix ⟨C in ⟩ and the perturbation term δC in (Eq. S22). In the following, we express theoretically the bias induced by this perturbation term on the estimation of T in . In particular, we will show how it scales with the parameter N W and the focusing quality. We consider here the input correlation guide star synthesized by a coherent combination of focal spots. By alternating aberration compensation at input and output, we can improve gradually the coherence factor C out/in and address forward multiple scattering associated with smaller isoplanatic patches (decrease N W ) while maintaining the bias δT in/out at a sufficiently low level. Figure S7 illustrates the importance of an alternate compensation of aberration and scattering at input and output. Figure S7b shows the evolution of the RPSF at each step of the algorithm when balancing between input and output. Figure S7c shows the evolution of the RPSF when the algorithm is only iterated at input. While a continuous balance between input and output aberration phase laws allows us to reach a diffraction-limited resolution at the end of the process (Fig. R3b), the absence of correction at output prevents from a refinement of the virtual guide star and does not allow our algorithm to converge towards a satisfying estimation of the matrix T in .

S12. CONVERGENCE OF THE MULTI-SCALE ANALYSIS PROCESS

The multi-scale process shown in Fig. S5 shall be stopped at some iterative step. Indeed, the spatial window W L cannot be reduced to a speckle grain otherwise the method would lead to a bucket image that consists in an incoherent summation of each de-scanned wave-field.

Qualitatively, the end of the process can be determined by a careful look at the image. An incoherent compensation of aberrations induces a loss of contrast on the final image. Figure S8 illustrates this assertion by comparing the original image (Fig. S8a), the RMI image obtained with a T-matrix of optimal resolution (6 × 6 µm 2 , see Fig. S8b) and a RMI image based on too small spatial windows W L (3 × 3 µm 2 , see Fig. S8c). The contrast of each image

] / ⟨I(ρ, z)⟩, tends to gradually increase when the estimator T approaches T (see comparison between Figs. S8a andb) and decrease when the compensation of aberrations and scattering becomes bucket-like (see comparison between Figs. S8b andc).

For the images displayed in Figs. S8a, b andc, we find F ∼ 1.48, F ∼ 1.61 and F ∼ 1.37, respectively. Nevertheless, an optimization criterion only based on the image contrast can be misleading since the contrast also depends on the sample reflectivity distribution.

A more reliable observable is the spatial correlation function C S (r p , r ′ p ) of the scattering component of the T-matrix between neighboring points r p and r ′ p (Methods). Examples of this spatial correlation function are displayed in Fig. S8d and Fig. S8e. While a spatial window of 6 × 6 µm 2 preserves a short-range correlation between neighbor windows (Fig. S8d), a spatial window of 3 × 3 µm 2 leads to a fully spatially incoherent estimator T (Fig. S8e). This observable clearly shows whether the estimator T leads to a coherent (i.e physical) or incoherent (i.e bucket-like) compensation of scattering. The number of iterations in the phase reversal algorithm has thus been based on this T-matrix correlation criterion. in contrast is explained by the fact that a part of the multiple scattering background is not addressed by RMI.

S14. DISCREPANCY BETWEEN INPUT AND OUTPUT T-MATRICES

FIG. S10. Discrepancy between input and output T-matrices. a. Scalar product P (ρ p , z)

between Tin and Tout at depth z =200 µm. b. Phase of the first pupil singular vector U 1 of Tin .

c. Scalar product P ′ (ρ p , z) between T′ in and Tout

While spatial reciprocity implies a strict equality between the wave distortions undergone by the incident and reflected waves in the sample arm of our experimental set-up, the input and output estimators of the T-matrix are far from it (see Figs. 4e andf). This discrepancy can be quantified by computing the normalized scalar product P (ρ p , z) between the coefficients of Tin and Tout :

Figure S10 shows the transverse evolution of this scalar product at depth z = 200 µm. As it could be anticipated when looking at the T-matrices in Figs. 4e andf, this scalar product is quite low: P (ρ p , z) = 0.3 in average. As we will see, this discrepancy can be, at least partially, explained by the aberrations in the reference arm. Indeed, according to Eq. S10, the input transmission matrix accumulates the aberrations undergone by the incident wave-field in the sample arm and the aberrations undergone by the reference wave-field,

On the contrary, Tout only grasps the wave distortions undergone by the reflected wave-field in the sample arm. If, in a first approximation, we assume that the aberration due to the

The SVD of T writes:

where σ k are the singular values arranged in decreasing order.

are the singular vectors of T in the pupil and focal plane, respectively. For a physical interpretation of these vectors, we take advantage of the equivalence between the SVD of T and the eigenvalue decomposition of the spatial correlation matrix,

The elements of C T correspond to the correlation coefficients between aberration phase laws obtained for each image pixel ρ p and ρ ′ p :

The first eigenvector V 1 of C T is thus the spatial domain where the degree of correlation between aberration phase laws is maximized. This degree of correlation is quantified by the normalized eigenvalue σ 2 1 , such that

The corresponding singular vector

is the transmittance the most-spatially invariant across the field-of-view. The same process can be iterated on the matrix T -σ 1 U 1 × V † 1 to retrieve the second eigenstate and so on. A set of orthogonal isoplanatic modes V k is finally obtained with a degree of correlation σ 2 k that decreases with their rank.

Figure S11 shows the result of the SVD of Tout for z = 50 µm (same depth as the one considered on top of Fig. 6 of the accompanying paper). Figure S11a displays its normalized singular values. A few predominant eigenvalues associated with the main isoplanatic modes seem to emerge from a continuum of lower eigenvalues associated with a multiple scattering background in each case. Figure S11c shows the four first eigenstates of Tout . While the first eigenstate V 1 spans over the whole field-of-view, the higher order isoplanatic modes V k are associated with spatial domains whose size decrease with the rank k of the eigenstate.

Another advantage of RMI versus AO consists in our ability of simulating any physical experiment in post-processing. If performed experimentally with an adaptive optics set up, the multi-scale compensation of wave distortions described in this paper would require:

(i) a complex adaptive optics arrangement to compensate for wave distortions both in the sample and reference arms; (ii) an extremely long acquisition time since the focusing process would have to be repeated 12 times on each of the 10 8 points of the imaged volume. The performance of matrix imaging is therefore impossible to reach with conventional adaptive optics tools.

S17. COMPARISON BETWEEN TIME REVERSAL AND ITERA-TIVE PHASE REVERSAL

Even though the ITR process does not provide a satisfying estimator at the first iteration of the matrix imaging process (see Supplementary Section S7 and Fig. S3), the multi-scale analysis and continuous balance between input and output correction enable to gradually reduce the size of the virtual scatterer and improves the ITR estimator at depth z =200 µm. In Figs. S13a 2 and a 3 , we compare the confocal images provided by ITR and IPR for L = 25 µm. Qualitatively, a careful look at the images shows that the ITR correction leads to a smoother image of the sample reflectivity than IPR. This is in agreement with Fig. S3 that showed that ITR tends to concentrate on a low spatial frequency spectrum while IPR addresses the whole spatial frequency components of the wave-field. In terms of image resolution, IPR thus shows a better performance than ITR.

In terms of contrast, the relative performance between IPR and ITR can be assessed by investigating the confocal energy at the end of the whole process. Fig. S13a 1 shows the ratio between their confocal gains at depth z =200 µm. The IPR process exhibits a better performance on a major part of the field-of-view. The ITR process is only better when it can hang on a highly reflecting structure. In that case, the maximization of backscattered energy on which the ITR process is based can lead to a better result than a criterion based on coherence as done by IPR. In practice, one could apply both IPR and ITR and keep the best option. Nevertheless, for sake of clarity and image continuity, the IPR process has been considered for all figures shown in the accompanying paper. Indeed, for small spatial windows (L = 6 µm), the ITR process leads to an image with strong vignetting to estimate input and output aberration phase laws. The CLASS algorithm amounts to find the aberration phase laws that maximize the confocal intensity in the focused basis. The input and output aberration laws are computed simultaneously and the whole process is then iterated several times to converge towards satisfying phase laws. On the contrary, our algorithm maximizes sequentially the coherence of input and output wave-fields generated by the virtual guide star. It allows us to decrease gradually the size of this guide star and improve the estimation of the T-matrix, while improving its resolution by gradually reducing the size of spatial windows. To compare the performance of both approaches, we have applied CLASS on our experimental data at z=200 µm. Figure S14 shows the results obtained by CLASS and IPR when the extension L of spatial windows W L is sufficiently large to allow the convergence of each method (L = 14 µm). Although the corrected RPSF display similar properties in average for each method [Figs. S14b 2 and c 2 ], the obtained images show significant differences [Figs. S14b 1 andc 1 ]. On the one hand, the CLASS image seems more contrasted because it tends to focus on the main scatterers of the field-of-view. On the other hand, the IPR