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 $\delta T_{in/out}$ at a sufficiently low level.

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

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}$.
<figure><img src="image_14.png" /><figcaption>FIG. S8. On the importance of a multi-scale compensation of wave distortions. a. Original confocal image of the cornea at *z* = 200 µm. b. Full-field correction. c. Direct compensation of aberrations over reduced spatial windows (*L* = 6 µm). d. Multi-scale strategy. Scale bar: 50 µm</figcaption></figure>

H. Validation of the method with a ground-truth object

We now provide an experimental validation of the method using a resolution target behind
the scattering medium. Although such a specular object does not reproduce the reflectivity
properties of tissues, this reference experiment will allow us to validate our multi-scale
analysis of wave distortions and also outline its limits.

This experiment is displayed in Supplementary Fig. S9. It consists in the imaging of a
resolution target placed right behind a 500-µm-thick mouse peritoneum layer (Supplementary
Fig. S9a). This tissue layer roughly displays the same scattering properties as the cornea
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², see Fig. S8b) and a RMI image based on too small spatial windows $W_L$ (3 × 3 µm², see Fig. S8c). The contrast of each image $I(\rho, z)$, $\mathcal{F}(z) = \text{std}[I(\rho, z)] / \langle I(\rho, z) \rangle$, tends to gradually increase when the estimator $\hat{\mathbf{T}}$ approaches **T** (see comparison between Figs. S8a and b) and decrease when the compensation of aberrations and scattering becomes bucket-like (see comparison between Figs. S8b and c). For the images displayed in Figs. S8a, b and c, we find $\mathcal{F} \sim 1.48$, $\mathcal{F} \sim 1.61$ and $\mathcal{F} \sim 1.37$, respectively. Nevertheless, an optimization criterion only based on the image contrast can be misleading since the contrast also depends on the sample reflectivity distribution.

A more reliable observable is the spatial correlation function $C_S(\mathbf{r}_p, \mathbf{r}_p')$ of the scattering component of the **T**-matrix between neighboring points $\mathbf{r}_p$ and $\mathbf{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² preserves a short-range correlation between neighbor windows (Fig. S8d), a spatial window of 3 × 3 µm² leads to a fully spatially incoherent estimator $\hat{\mathbf{T}}$ (Fig. S8e). This observable clearly shows whether the estimator $\hat{\mathbf{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.
used in the accompanying manuscript.

<figure><img src="image_15.png" /><figcaption>FIG. S9. **Imaging a resolution target hidden behind an opaque tissue layer.** **a.** Experimental configuration. Image credit: Setreset on Wikimedia Commons. **b.** Incoherent image (colorbar in linear scale). **c.** FFOCT image (B&W bar in log-scale). **d.** Matrix Image (B&W bar in log-scale). **e-h.** Zooms on the smallest details of the resolution target images displayed in panels **b-d**, respectively. Scale bar: 50 µm.</figcaption></figure>

Supplementary Figure S9b displays the incoherent image of the target in reflection which is obtained here by blocking the reference arm in the experimental set up. Its foggy feature highlights the strong turbidity of the scattering layer. Supplementary Figure S9c shows the FFOCT image acquired for $\Delta\rho_{in} = 0$. Its comparison with Supplementary Fig. S9b illustrates the drastic filtering of diffuse multiple scattering operated by the time gating process in FFOCT. It also shows how the confocal filter allows to reveal the large patterns of the resolution target. The fact that FFOCT is robust with respect to aberrations for specular objects has already been noticed in a previous work$^{29}$. In that case, aberrations only give rise to fluctuations of image contrast across the field-of-view. Nevertheless, this robustness to aberrations vanishes for the smallest details of the target that cannot be detected [see the corresponding zoom in Supplementary Fig. S9g].

The aberration and scattering induced by the tissue layer can be evaluated by the RPSF whose spatial evolution is displayed in Supplementary Fig. S10a. It shows drastic variations across the field-of-view, a manifestation of a particularly short-range memory effect across
<figure><img src="image_15.png" /><figcaption>FIG. S8. **Confocal images at several steps of the multi-scale analysis.** a Initial en-face image of the cornea at depth $z \sim 100$ µm. b,c RMI images based on a T-matrix estimator of spatial resolution $L = 6$ µm and $L = 3$ µm, respectively. d,e Spatial correlation $C_S$ of $\hat{T}^S$ with respect to one reference location (white arrow). Scale bars : 50 µm.</figcaption></figure>
the field-of-view. Compared with the incoherent RPSFs measured inside the cornea (Fig. 4c of the accompanying paper), the coherent RPSFs here show a more contrasted feature. This is due to the specular nature of the sample. In that regime, the RPSF is a coherent function of the focusing quality since it scales as the convolution between the coherent input and output PSFs³³:

$$I(\Delta\rho_{\text{in}}, \rho_{\text{out}}, z) = |H_{\text{in}}^{(l)} \overset{\Delta\rho_{\text{in}}}{\otimes} H_{\text{out}}^{(l)}(\Delta\rho_{\text{in}}, \rho_{\text{out}}, z)|^2.$$

On the contrary, for a sample of random reflectivity (like the opaque cornea under study), the RPSF is an incoherent measure of the focusing quality. Its ensemble average can be expressed as the convolution between the incoherent input and output PSFs (Eq. 3 of the accompanying paper):

$$\langle I(\Delta\rho_{\text{in}}, \rho_{\text{out}}, z) \rangle \propto |H_{\text{in}}^{(l)}|^2 \overset{\Delta\rho_{\text{in}}}{\otimes} |H_{\text{out}}^{(l)}|^2 (\Delta\rho_{\text{in}}, \rho_{\text{out}}, z).$$

Despite the short-range memory effect highlighted by the original RPSFs (Supplementary Fig. S10a), a multi-scale compensation of aberration and scattering phenomena allows us to retrieve an almost diffraction-limited RPSF across the whole surface (Supplementary Fig. S10b), except inside the patterns of the resolution target since there is no back-scattered wave-field there. The resulting image is displayed in Supplementary Fig. S9d. Compared with its original counterpart (Supplementary Fig. S9c), an homogeneous contrast is obtained throughout the field-of-view. Above all, RMI is able to retrieve the smallest details of the resolution target (Supplementary Fig. S9g) that FFOCT initially failed to reveal (Supplementary Fig. S9f). Nevertheless, the image is not perfect. A first reason comes from the limited de-scan of the RPSF. This loss of information gives rise to a residual incoherent background in the final RPSF (Supplementary Fig. S10c). The second reason is the limited isoplanicity. The size L of spatial windows at the end of the RMI process is limited to 10 µm. Hence we cannot compensate for scattering phenomena giving rise to a memory effect range smaller than L.
S13. QUANTIFYING THE CONTRAST ENHANCEMENT

<figure><img src="image_16.png" /><figcaption>FIG. S9. **Confocal gain provided by the matrix imaging process.** ab. Transverse cross-section 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.</figcaption></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 and f) 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 and f, we could have expected a larger confocal intensity enhancement. This moderate gain
<figure><img src="image_16.png" /><figcaption>FIG. S10. Reflection point spread function in the resolution target experiment. a. Map of initial RPSFs. b. Map of final RPSFs. c. Mean RPSF after aberration compensation (scale bar: 2 µm).</figcaption></figure>
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

<figure><img src="image_17.png" /><figcaption>FIG. S10. Discrepancy between input and output T-matrices. a. Scalar product $P(\rho_p, z)$ between $\hat{T}_{in}$ and $\hat{T}_{out}$ at depth $z = 200$ µm. b. Phase of the first pupil singular vector $U_1$ of $\hat{T}_{in}$. c. Scalar product $P'(\rho_p, z)$ between $\hat{T}'_{in}$ and $\hat{T}_{out}$</figcaption></figure>

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 and f). This discrepancy can be quantified by computing the normalized scalar product $P(\rho_p, z)$ between the coefficients of $\hat{T}_{in}$ and $\hat{T}_{out}$:

$$P(\rho_p, z) = N_u^{-1} \sum_{\mathbf{u}} \hat{T}_{\text{in}}(\mathbf{u}, \rho_p, z) \hat{T}_{\text{out}}^*(\mathbf{u}, \rho_p, z) \quad (\text{S50})$$

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 and f, this scalar product is quite low: $P(\rho_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,

$$\mathbf{T}_{\text{in}} = \mathbf{T} \circ \mathbf{T}_{\text{ref.}} \quad (\text{S51})$$

On the contrary, $\hat{T}_{out}$ 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
## S4. ANALYSIS OF THE TRANSMITTANCE MATRIX

### A. Discrepancy between input and output $T$-matrices

<figure><img src="image_17.png" /><figcaption>FIG. S11. Discrepancy between input and output $T$-matrices. a. Scalar product $P(\rho_p, z)$ between $\hat{T}_{in}$ and $\hat{T}_{out}$ at depth $z=200$ µm. b. Phase of the first pupil singular vector $U_1$ of $\hat{T}_{in}$. c. Scalar product $P'(\rho_p, z)$ between $\hat{T}'_{in}$ and $\hat{T}_{out}$. Scale bar: 50 µm.</figcaption></figure>

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 and f of the accompanying paper). This discrepancy can be quantified by computing the normalized scalar product $P(\rho_p, z)$ between the coefficients of $\hat{T}_{in}$ and $\hat{T}_{out}$:

$$P(\rho_p, z) = N_u^{-1} \sum_{\mathbf{u}} \hat{T}_{in}(\mathbf{u}, \rho_p, z) \hat{T}_{out}^*(\mathbf{u}, \rho_p, z) \quad (S48)$$

Supplementary Fig. S11 shows the transverse evolution of this scalar product at depth $z = 200$ µm. As it could be anticipated when looking at the $\hat{T}$-matrices in Figs. 4e and f, this scalar product is quite low: $P(\rho_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. S8, 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,

$$T_{in} = T \circ T_{ref}. \quad (S49)$$

On the contrary, $\hat{T}_{out}$ 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
reference arm is isoplanatic, it can be extracted by considering the first eigenstate of $\hat{T}_{\text{in}}$ (see Supplementary Section S15). The phase of the pupil singular vector $\mathbf{U}_1$ displayed in Fig. S10 is an estimator of $\mathbf{T}_{\text{ref}}$. Not surprisingly, it mainly corresponds to a spherical aberration phase law. One can subtract this reference phase to $\hat{T}_{\text{in}}$ in order to build a matrix $\hat{T}'_{\text{in}} = \hat{T}_{\text{in}} \circ \mathbf{U}_1^*$. One can expect the scalar product $P'(\rho_p, z)$ between $\hat{T}'_{\text{in}}$ and $\hat{T}_{\text{out}}$ to be increased compared to its initial value $P(\rho_p, z)$ (see comparison between Figs. S10a and c). 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 $\hat{T}_{\text{in}}$ and $\hat{T}_{\text{out}}$.

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³⁸, 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.
reference arm is isoplanatic, it can be extracted by considering the first eigenstate of $\hat{T}_{\text{in}}$
(see Supplementary Section S4 B). The phase of the pupil singular vector $\mathbf{U}_1$ displayed in
Supplementary Fig. S11 is an estimator of $\mathcal{T}_{\text{ref}}$. Not surprisingly, it mainly corresponds to
a spherical aberration phase law. One can subtract this reference phase to $\hat{T}_{\text{in}}$ in order to
build a matrix $\hat{T}'_{\text{in}} = \hat{T}_{\text{in}} \circ \mathbf{U}_1^*$. One can expect the scalar product $P'(\rho_p, z)$ between $\hat{T}'_{\text{in}}$
and $\hat{T}_{\text{out}}$ to be increased compared to its initial value $P(\rho_p, z)$ (see comparison between
Supplementary Figs. S11a and c). This is actually what we observe even though the scalar
product $P'$ remains smaller than 0.7 (Supplementary Fig. S11c). It means that the spherical
aberration law induced by the reference arm account partially for the mismatch between $\hat{T}_{\text{in}}$
and $\hat{T}_{\text{out}}$.

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 S3 C.

B. 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<sup>40</sup>, 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.
The SVD of $\hat{T}$ writes:

$$ \hat{T} = \sum_k \sigma_k \mathbf{U}_k \mathbf{V}_k^\dagger \qquad (S52) $$

where $\sigma_k$ are the singular values arranged in decreasing order. $\mathbf{U}_k = [U_k(\mathbf{u})]$ and $\mathbf{V}_k = [V_k(\boldsymbol{\rho}_p)]$ are the singular vectors of $\hat{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 $\hat{T}$ and the eigenvalue decomposition of the spatial correlation matrix,

$$ \mathbf{C}_T = \hat{T} \hat{T}^\dagger. \qquad (S53) $$

The elements of $\mathbf{C}_T$ correspond to the correlation coefficients between aberration phase laws obtained for each image pixel $\boldsymbol{\rho}_p$ and $\boldsymbol{\rho}_p'$:

$$ C_T(\boldsymbol{\rho}_p, \boldsymbol{\rho}_p') = \sum_{\mathbf{u}} \hat{T}(\mathbf{u}, \boldsymbol{\rho}_p) \hat{T}^*(\mathbf{u}, \boldsymbol{\rho}_p). \qquad (S54) $$

The first eigenvector $\mathbf{V}_1$ of $\mathbf{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 $\bar{\sigma}_1^2$, such that

$$ \bar{\sigma}_k^2 = \frac{\sigma_k^2}{\sum_l \sigma_l^2} = \frac{\mathbf{V}_k^\dagger \times \mathbf{C}_T \times \mathbf{V}_k}{\mathrm{Tr}\{\mathbf{C}_T\}} \qquad (S55) $$

The corresponding singular vector

$$ \mathbf{U}_1 = \sigma_1^{-1} \hat{T} \times \mathbf{V}_1 \qquad (S56) $$

is the transmittance the most-spatially invariant across the field-of-view. The same process can be iterated on the matrix $\hat{T} - \sigma_1 \mathbf{U}_1 \times \mathbf{V}_1^\dagger$ to retrieve the second eigenstate and so on. A set of orthogonal isoplanatic modes $\mathbf{V}_k$ is finally obtained with a degree of correlation $\bar{\sigma}_k^2$ that decreases with their rank.

Figure S11 shows the result of the SVD of $\hat{T}_{\text{out}}$ 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 $\hat{T}_{\text{out}}$. While the first eigenstate $\mathbf{V}_1$ spans over the whole field-of-view, the higher order isoplanatic modes $\mathbf{V}_k$ are associated with spatial domains whose size decrease with the rank $k$ of the eigenstate.
The SVD of $\hat{\mathcal{T}}$ writes:

$$ \hat{\mathcal{T}} = \sum_k \sigma_k \mathbf{U}_k \mathbf{V}_k^\dagger \qquad (S50) $$

where $\sigma_k$ are the singular values arranged in decreasing order. $\mathbf{U}_k = [U_k(\mathbf{u})]$ and $\mathbf{V}_k = [V_k(\boldsymbol{\rho}_p)]$ are the singular vectors of $\hat{\mathcal{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 $\hat{\mathcal{T}}$ and the eigenvalue decomposition of the spatial correlation matrix,

$$ \mathbf{C}_{\mathcal{T}} = \hat{\mathcal{T}} \times \hat{\mathcal{T}}^{\dagger}. \qquad (S51) $$

The elements of $\mathbf{C}_{\mathcal{T}}$ correspond to the correlation coefficients between aberration phase laws obtained for each image pixel $\boldsymbol{\rho}_p$ and $\boldsymbol{\rho}_p'$:

$$ C_{\mathcal{T}}(\boldsymbol{\rho}_p, \boldsymbol{\rho}_p') = \sum_{\mathbf{u}} \hat{\mathcal{T}}(\mathbf{u}, \boldsymbol{\rho}_p) \hat{\mathcal{T}}^*(\mathbf{u}, \boldsymbol{\rho}_p). \qquad (S52) $$

The first eigenvector $\mathbf{V}_1$ of $\mathbf{C}_{\mathcal{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 $\bar{\sigma}_1^2$, such that

$$ \bar{\sigma}_k^2 = \frac{\sigma_k^2}{\sum_l \sigma_l^2} = \frac{\mathbf{V}_k^\dagger \times \mathbf{C}_{\mathcal{T}} \times \mathbf{V}_k}{\mathrm{Tr}\{\mathbf{C}_{\mathcal{T}}\}} \qquad (S53) $$

The corresponding singular vector

$$ \mathbf{U}_1 = \sigma_1^{-1} \hat{\mathcal{T}} \times \mathbf{V}_1 \qquad (S54) $$

is the transmittance the most-spatially invariant across the field-of-view. The same process can be iterated on the matrix $\hat{\mathcal{T}} - \sigma_1 \mathbf{U}_1 \times \mathbf{V}_1^\dagger$ to retrieve the second eigenstate and so on. A set of orthogonal isoplanatic modes $\mathbf{V}_k$ is finally obtained with a degree of correlation $\bar{\sigma}_k^2$ that decreases with their rank.

Supplementary Figure S12 shows the result of the SVD of $\hat{\mathcal{T}}_{\mathrm{out}}$ for $z = 50$ µm (same depth as the one considered on top of Fig.6 of the accompanying paper). Supplementary Figure S12a 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. Supplementary Figure S12c shows the four first eigenstates of $\hat{\mathcal{T}}_{\mathrm{out}}$. While the first eigenstate $\mathbf{V}_1$ spans over the whole field-of-view, the higher order isoplanatic modes $\mathbf{V}_k$ are associated with spatial domains whose size decrease
<figure><img src="image_18.png" /><figcaption>FIG. S11. Singular value decomposition of the T-matrix at depth $z = 50$ µm. a Normalized singular values of $\hat{T}_{\text{out}}$. b Degree of vorticity associated with each pupil singular vector $U_k$ of $\hat{T}_{\text{out}}$. c Four first eigenstates of $\hat{T}_{\text{out}}$: Phase of the transmittance $U_k$ and modulus of the corresponding isoplanatic modes $V_k$.</figcaption></figure>

The complexity (i.e the spatial frequency content) of the associated transmittance $U_k$ also increases with this rank.

The nature of the associated wave distortions can be investigated by considering the phase $\phi_k$ of each singular vector $U_k$. More precisely, recent works$^{41,42}$ showed how aberrations and scattering can be discriminated by computing the divergence and curl of the phase gradient $\nabla\phi_k$. Each phase law $\phi_k$ can be decomposed into: (i) an irrotational component $\phi_k^{(\text{irr})}$, such that $\nabla \wedge \nabla\phi_k^{(\text{irr})} = 0$, associated with low-order aberrations; (ii) a curl component
<figure><img src="image_18.png" /><figcaption>FIG. S12. Singular value decomposition of the $\hat{T}$-matrix at depth $z = 50$ µm. a. Normalized singular values of $\hat{T}_{\text{out}}$. b. Degree of vorticity associated with each pupil singular vector $U_k$ of $\hat{T}_{\text{out}}$. c. Four first eigenstates of $\hat{T}_{\text{out}}$: Phase of the transmittance $U_k$ and modulus of the corresponding isoplanatic modes $V_k$.</figcaption></figure>

with the rank $k$ of the eigenstate. The complexity (i.e the spatial frequency content) of the associated transmittance $U_k$ also increases with this rank.

The nature of the associated wave distortions can be investigated by considering the phase $\phi_k$ of each singular vector $U_k$. More precisely, recent works$^{42,60}$ showed how aberrations and scattering can be discriminated by computing the divergence and curl of the phase gradient $\nabla\phi_k$. Each phase law $\phi_k$ can be decomposed into: (i) an irrotational component $\phi_k^{(\text{irr})}$, such that $\nabla \wedge \nabla\phi_k^{(\text{irr})} = 0$, associated with low-order aberrations; (ii) a curl component
$\phi_k^{(\text{rot})}$, such that $\nabla \cdot \nabla \phi_k^{(\text{rot})} = 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 $\alpha$ of vorticity can be assessed by looking at the ratio between the energy of each contribution, such that:

$$ \alpha_k = \frac{|\nabla \wedge \nabla \phi_k|^2}{|\Delta \phi_k|^2} \qquad (S57) $$

This degree of vorticity is displayed for each eigenstate in Fig. S11b. Although it shows some fluctuations, $\alpha_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.

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$) 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$).
$\phi_k^{(\text{rot})}$, such that $\nabla \cdot \nabla \phi_k^{(\text{rot})} = 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 $\alpha$ of vorticity can be assessed by looking at the ratio between the energy of each contribution, such that:

$$ \alpha_k = \frac{|\nabla \wedge \nabla \phi_k|^2}{|\Delta \phi_k|^2} \qquad (S55) $$

This degree of vorticity is displayed for each eigenstate in Supplementary Fig. S12b. Although it shows some fluctuations, $\alpha_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 can be a manifestation of forward multiple scattering paths.

### C. Angular decomposition of the $\mathcal{T}$-matrix

To investigate the effect of forward multiple scattering, an angular decomposition of the transmitted wave-field between the cornea surface and the focal plane can be performed. Interestingly, this can be done by considering the Fourier transform of the transmission matrix estimator $\tilde{\mathcal{T}}(z) = [\hat{\mathcal{T}}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z)]$ along the focused basis:

$$ \tilde{\mathcal{T}}(z) = \hat{\mathcal{T}}(z) \times \mathbf{T}_0, \qquad (S56) $$

which writes, in terms of matrix coefficients,

$$ \tilde{\mathcal{T}}(\mathbf{u}, \Delta\mathbf{u}, z) = \sum_{\boldsymbol{\rho}_{\text{out}}} \hat{\mathcal{T}}](\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) \exp \left( -i \frac{2\pi}{\lambda f} \Delta\mathbf{u} \cdot \boldsymbol{\rho}_{\text{out}} \right). \qquad (S57) $$

To show the relationship between $\tilde{\mathcal{T}}(\mathbf{u}_{\text{in}}, \Delta\mathbf{u}, z)$ and the angular distribution of light in the focal plane, one can use the relationship between the transmission and transmittance matrices,

$$ \mathbf{T} = \mathcal{T} \circ \mathbf{T}_0, \qquad (S58) $$

which writes in terms of matrix coefficients:

$$ T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = \mathcal{T}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) \exp \left( j \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_{\text{out}} \right), \qquad (S59) $$

or, equivalently,

$$ \mathcal{T}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) \exp \left( -j \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_{\text{out}} \right). \qquad (S60) $$
<figure><img src="image_19.png" /><figcaption>FIG. S12. Decomposition of the focusing law $\phi$ 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 bar; 50 µm].

**c** Associated confocal gain. **c** Example of local RPSF obtained in the white rectangle displayed in panel **b**1).</figcaption></figure>
Injecting this last expression into Eq. S58 leads to the following expression for $\tilde{T}$:

$$ \tilde{T}(\mathbf{u}_{\text{in}}, \Delta\mathbf{u}, z) = \sum_{\Delta\mathbf{u}} \exp \left[ -j \frac{2\pi}{\lambda f} (\Delta\mathbf{u} + \mathbf{u}_{\text{in}}) \cdot \boldsymbol{\rho}_{\text{out}} \right] T(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}) \quad (\text{S61}) $$

$$ = \tilde{T}(\mathbf{u}_{\text{in}}, \mathbf{u}_{\text{in}} + \Delta\mathbf{u}). \quad (\text{S62}) $$

Hence, the $\tilde{T}$-matrix actually corresponds to the de-scan transmission matrix in the pupil plane. Its coefficients $\tilde{T}(\mathbf{u}_{\text{in}}, \Delta\mathbf{u}, z)$ provide the angular dispersion of the transmitted wave-field at depth $z$ with respect to the incident plane wave of transverse wave vector $\mathbf{k}_{\text{in}} = \mathbf{u}_{\text{in}}/f$ (Supplementary Fig. S13a). For small angles, the deviation angle $\Delta\theta = (\Delta\theta_x, \Delta\theta_y)$ can be expressed as follows:

$$ \Delta\theta \sim \sin\theta_{\text{out}} - \sin\theta_{\text{in}} = \Delta\mathbf{u}/f \quad (\text{S63}) $$

The angular dispersion of the wave-field between the cornea surface and the focal plane can thus be obtained by averaging the intensity of the transmitted wave-field over the incident wave vector:

$$ P_{\mathcal{T}}(\Delta\mathbf{u}, z) = \langle |\tilde{T}(\mathbf{u}_{\text{in}}, \Delta\mathbf{u}/f, z)|^2 \rangle_{\mathbf{u}_{\text{in}}}. \quad (\text{S64}) $$

Supplementary Figure S13b shows this angular distribution at depths $z = 50$ and $200$ µm, respectively. It displays the following shape: A close-to-ballistic peak around $\Delta\mathbf{u} = 0$ on top of a wide pedestal generated by forward multiple scattering. This statement is confirmed by investigating the angular distributions, $P_A(\Delta\mathbf{u}, z)$ (Supplementary Fig. S13c) and $P_S(\Delta\mathbf{u}, z)$ (Supplementary Fig. S13d), associated with the aberration and scattering matrices, $\mathcal{A}_{\text{out}}$ and $\mathcal{S}_{\text{out}}$, respectively. As anticipated, the aberration component of the $\mathcal{T}$-matrix is associated with a close-to-ballistic peak around $\Delta\mathbf{u} = 0$ (Supplementary Fig. S13c), while its scattering component gives rise to wide distribution of deviation angles in the focal plane (Supplementary Fig. S13d). As illustrated by Supplementary Fig. S13h, a wider angular distribution is observed for $P_S$ at $200$ µm than at $50$ µm. The angular width of the photon distribution at -10 dB goes from $4\text{deg}$ at $z = 50\mu\text{m}$ to $10\text{deg}$ at $z = 200\mu\text{m}$. This angular dispersion of the transmitted wave-field between the cornea surface and the focal plane is a manifestation of scattering events taking place between those two planes. It thus confirms that our $\mathcal{T}$-matrix estimator indeed contains a forward multiple scattering contribution.
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 ITERATIVE TIME REVERSAL AND ITERATIVE 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₂ and a₃, 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₁ 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
<figure><img src="image_19.png" /><figcaption>FIG. S13. Plane wave decomposition of the **T**-matrix. **a.** Definition of the incident and transmitted wave vectors, **k**<sub>in</sub> and **k**<sub>out</sub>, and their associated angles, θ<sub>in</sub> and θ<sub>out</sub>, at the cornea surface and the focal plane. **b-d.** Wave vector deviation distribution of the transmitted wave-field P<sub>T</sub>(Δ**u**, z) (**b**), of its aberrated component P<sub>A</sub>(Δ**u**, z) (**c**), and of its scattered component P<sub>S</sub>(Δ**u**, z) (**d**) [Top: z = 50 µm; Bottom: z = 200 µm]. **e.** Angular distribution P<sub>S</sub>(Δθ, z) of the scattered wave-field at depths z= 50 µm (blue line) and z= 200 µm (red line).</figcaption></figure>

**D. Coherent backscattering as a manifestation of multiple scattering**

Another key observable to prove the existence of multiple scattering paths is the coherent back-scattering phenomenon. This phenomenon results from the constructive interference between reciprocal multiple scattering paths. Investigated originally in a plane wave basis (k-space)<sup>34,35,61</sup>, it manifests as an enhancement by a factor two of the time-gated intensity in the back-scattering direction<sup>62</sup>. Studied in the focused basis (position space), it gives rise to an enhancement by a factor two of the mean intensity at Δρ = 0<sup>33,63,64</sup> (Supplementary Fig. S14). This phenomenon is investigated in Supplementary Fig. S15. The initial maps of RPSFs provided by Fig. 4c and reproduced in Supplementary Fig. S15a does not exhibit a clear signature of the CBS because of the aberrations induced by our imaging system that alters its shape<sup>65</sup>.

After compensation of the reference arm aberration (Supplementary Fig. S11b) at the
first step of our aberration correction process, a new map of RPSFs is obtained and shown
<figure><img src="image_20.png" /><figcaption>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.</figcaption></figure>

effects (Fig. S13b<sub>3</sub>), which is a manifestation of a lack of correlations of the estimator $\hat{T}$ between adjacent windows. On the contrary, the IPR leads to an estimator whose scattered component $\hat{S}$ 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<sup>20,23,24</sup>. Based on the recording of a time-gated reflection matrix, it exploits its input-output correlations in the plane wave basis
<figure><img src="image_20.png" /></figure>

FIG. S14. **Coherent back-scattering phenomenon in the focused basis.** The RPSF is made of two contributions in the multiple scattering regime: **a** A diffuse component which results from the self-interference of each scattering path with itself; **b** A coherent intensity resulting from the interference between reciprocal multiple scattering paths inside the medium. Those paths undergo the same scattering sequence but in reverse order. **c** When the input and output focused beams coincide ($ρ_{in} = ρ_{out}$), the interference is constructive and leads to an enhancement by a factor two of the RPSF with respect to the diffuse background. This is the so-called coherent backscattering peak highlighted by Supplementary Fig. S15c.

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

FIG. S15. **Reflection point spread function in the focused basis.** **a.** Maps of initial RPSFs at $z = 200$ µm. **b.** Maps of RPSFs at the same depth after compensation of the isoplanatic aberration phase law (Fig. S10b) mainly induced by the reference arm. **c.** Radial average of the RPSF averaged over the area surrounded by a black square in panel **b**.

in Fig. S15b. It exhibits a following shape: A confocal peak due to the single scattering contribution and to the coherent back-scattering phenomenon that results from the constructive interference between reciprocal multiple scattering paths inside the medium (Supplementary Fig. S14c); on top of a diffusive background resulting from the incoherent summation of each
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.

<figure><img src="image_21.png" /><figcaption>FIG. S14. Comparison between CLASS and RMI upon convergence. a. Original confocal image ($\mathbf{a}_1$) and mean RPSF ($\mathbf{a}_2$) at $z=200$ µm. b, c. Corresponding CLASS ($\mathbf{b}_1$) and IPR ($\mathbf{c}_1$) images with their mean RPSFs ($\mathbf{b}_2$ and $\mathbf{c}_2$, respectively). The size of spatial windows is $L = 13$ µm.</figcaption></figure>

To compare the performance of both approaches, we have applied CLASS on our ex-
perimental 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 and c_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
multiple scattering path intensity (Supplementary Fig. S14a).

In absence of noise, the confocal intensity is therefore equal to $I_C = I_S + 2I_M$, with $I_S$ the single scattering intensity and $I_M$ the multiple scattering intensity. The incoherent background directly provides the multiple scattering intensity. The ratio $\beta_C = I_C/I_M$ between the confocal peak and the incoherent background can actually provide an estimate for the multiple scattering rate³³:

$$ \beta_M = I_M/(I_S + I_M) = 1/(\beta_C - 1) \qquad (S65) $$

or, equivalently, of the single scattering rate

$$ \beta_S = I_S/(I_S + I_M) = (\beta_C - 2)/(\beta_C - 1) \qquad (S66) $$

Note that, in practice, a quantitative measurement of the single/multiple scattering rates is not so easy to perform since the multiple scattering background does not exhibit a flat profile especially at shallow depth⁶⁶.

Anyway, a value of $\beta_C$ close to 2 means a predominant multiple scattering contribution. This is actually what we observe in many parts of the field-of-view, in particular in the area surrounded by a black rectangle in Supplementary Fig. S15b. The radial distribution of the RPSF displayed in Supplementary Fig. S15c shows a CBS enhancement of two, proof that multiple scattering is predominant in this region. At the end of the multi-scale aberration correction process, the maps of RPSF shows a much larger confocal ratio $\beta_C$, proof that multiple scattering trajectories have been (at least partially) compensated by our $\mathcal{T}$-matrix estimator.

## E. Measuring the Scattering Mean Free Path

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

Unfortunately, here, the cornea is not healthy but oedematous. The depth evolution of the confocal intensity in the stroma is thus strongly impacted by multiple scattering and cannot be used for a measurement of $\ell_s$. Moreover, in the epithelium, the different layers of
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.

<figure><img src="image_22.png" /><figcaption>FIG. S15. Comparison between CLASS and RMI for small isoplanatic patches. a. Original confocal image ($\mathbf{a}_1$) and mean RPSF ($\mathbf{a}_2$) at $z=200~\mu m$. b, c. Corresponding CLASS ($\mathbf{b}_1$) and IPR ($\mathbf{c}_1$) images with their mean RPSFs ($\mathbf{b}_2$ and $\mathbf{c}_2$, respectively). The size of spatial windows is $L = 6~\mu m$.</figcaption></figure>

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~\mu m$). A local RPSF obtained via CLASS is displayed in Fig. S15b₂. 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₂). When looking into details to the CLASS image (Fig. S15b₁), strong vignetting effects can be observed while the IPR image exhibits a continuous reflectivity (Fig. S15c₁).

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 $\mathbf{r}_p$ of the area
<figure><img src="image_22.png" /><figcaption>FIG. S16. Logarithm of the single scattering rate ln $\beta_S$ versus depth (blue dots) fitted with Eq. S67 (red dashed line).</figcaption></figure>

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<sup>69</sup>. It consists in investigating the depth evolution of the single scattering rate. In 3D, radiative transfer solution indicates the following exponential scaling for $\beta_S$<sup>36</sup>:

$$ \beta_S(z) \sim \exp(-z/\ell_s) \qquad (S67) $$

To estimate $\beta_S$ (Eq. S66), the confocal ratio $\beta_c$ has been estimated as follows:

$$ \hat{\beta}_C(z) = \frac{\max_{\Delta\rho}\{I(\Delta\rho, z)\}}{\min_{\Delta\rho}\{I(\Delta\rho, z)\}} \qquad (S68) $$

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

Supplementary Figure S16 displays the depth evolution of the single scattering rate $\beta_S(z)$ computed from $\hat{\beta}_C$ (Eq. S66). It exhibits an exponential decay in the stroma beyond $z = 140$ µm. The single scattering rate cannot be estimated beyond $z = 200$ µm because our estimator
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 $\hat{S}$ (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. S16 accounts for the vignetting effects observed for the CLASS image in Fig. S15b₁.

<figure><img src="image_23.png" /><figcaption>FIG. S16. Memory effect exhibited by the forward multiple scattering estimator: Comparison between CLASS and RMI. Spatial correlation function $C_S(\mathbf{r}_p, \mathbf{r}'_p)$ of the scattering component $\hat{S}$ for a given point $\mathbf{r}'_p$. **a.** CLASS. **b.** RMI. The depth is $z = 200$ µm and the size of spatial windows is $L = 6$ µm.</figcaption></figure>
of $\beta_C(z)$ starts to be impacted by the experimental noise. Therefore, the fit of $\beta_S(z)$ with Eq. S67 is performed from $z = 140$ to $z = 200$ µm. We find $\ell_s \sim 35$ µm.

### F. Quantifying the contrast enhancement

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

FIG. S17. Confocal gain provided by the matrix imaging process. a-b. Transverse cross-section of the confocal gain observed for the en-face images displayed in Fig. 3b at depths 50 µm, 250 µm and 350 µ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.

Supplementary Figure S17 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 and f) and that RMI tends to realign in phase by means of a digital optical phase conjugation. Supplementary Figure S17b clearly shows that the confocal gain increases with depth $z$. Indeed, multiple scattering becomes predominant in depth and the
S19. GLOSSARY

<table><caption>TABLE S1. Glossary of the operators used in this study</caption><thead><tr><th>operator</th><th>definition</th></tr></thead><tbody><tr><td><b>R</b> = [R(ρ<sub>in</sub>, ρ<sub>out</sub>, z)]</td><td>reflection matrix in the focused basis</td></tr><tr><td><b>R</b><sub>in/out</sub> = [R<sub>in</sub>(Δρ<sub>in/out</sub>, ρ<sub>out/in</sub>, z)]</td><td>de-scanned matrix at input/output</td></tr><tr><td><b>D</b><sub>in/out</sub> = [D(<b>u</b><sub>in/out</sub>, ρ<sub>out/in</sub>), z]</td><td>input/output distortion matrix</td></tr><tr><td><b>H</b><sub>in/out</sub> = [H<sub>in</sub>(ρ<sub>s</sub>, ρ<sub>in/out</sub>, z)]</td><td>input/output focusing matrix</td></tr><tr><td><b>F</b> = [F(<b>u</b>, Δρ)]</td><td>Fourier transform operator</td></tr><tr><td><b>C</b><sub>in/out</sub> = [C<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>u'</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output pupil correlation matrix of <b>D</b><sub>in</sub></td></tr><tr><td><b>C</b><sub>H</sub> = [C<sub>H</sub>(<b>u</b>, <b>u'</b>)]</td><td>correlation matrix associated with<br>a virtual reflector of reflectivity |H(ρ<sub>s</sub>)|^2</td></tr><tr><td><b>T</b><sub>in/out</sub> = [T<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output transmission matrix</td></tr><tr><td><b>T̂</b><sub>in/out</sub> = [T<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>estimator of the input/output transmission matrix</td></tr><tr><td>ϕ<sub>in/out</sub> = [ϕ<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>phase of the input/output transmission matrix estimator</td></tr><tr><td><b>A</b><sub>in/out</sub> = [T<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output aberration matrix</td></tr><tr><td><b>S</b><sub>in/out</sub> = [S<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output forward multiple scattering matrix</td></tr></tbody></table>
transmission phase laws become more and more complex. Note, however, that given the
complexity of phase laws displayed in Figs. 4e and f, we could have expected a larger confocal
intensity enhancement. This moderate gain in contrast is explained by the fact that a part
of the multiple scattering background is not addressed by RMI.

S5. COMPARISON WITH STATE-OF-THE-ART METHODS

**A. 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<sup>42</sup>. On the contrary,
the interferometric measurement of the reflected wave-field gives access to this scattering
component which is crucial for deep imaging.

Supplementary Figure S18 illustrates this assertion by first showing the decomposition
of the input and output phase laws (Supplementary Figs. S18a₁) into their irrotational
(Supplementary Figs. S18a₂) and curl (Supplementary Figs. S18a₃) 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 (Supplementary Fig. S18b). This
can be quantified by the confocal gain exhibited at the end of the matrix imaging process
(Supplementary Fig. S18c) and the corresponding RPSF (Fig. S18d). While the access to
the curl component of the focusing laws allows us to reach a confocal gain up to 13 dB
(Supplementary Fig. S18c₁), conventional AO would only allow a compensation of low-order
aberrations, giving rise to a weak confocal gain (< 3 dB, Supplementary Fig. S18c₂).

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
[1] B. Gérardin, J. Laurent, A. Derode, C. Prada, and A. Aubry, Phys. Rev. Lett. **113**, 173901 (2014).

[2] M. Horodynski, M. K'uhmayer, C. Ferise, S. Rotter, and M. Davy, Nature **607**, 281 (2022).

[3] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Nature **453**, 891 (2008).

[4] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nat. Phys. **4**, 945 (2008).

[5] V. Ntziachristos, Nat. Methods **7**, 603 (2010).

[6] J. Bertolotti and O. Katz, Nat. Phys. **18**, 1008 (2022).

[7] C. Dunsby and P. French, J. Phys. D: Appl. Phys. **36**, R207 (2003).

[8] A. Badon, A. C. Boccara, G. Lerosey, M. Fink, and A. Aubry, Opt. Express **25**, 28914 (2017).

[9] D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, C. Puliafito, and a. et, Science **254**, 1178 (1991).

[10] H. W. Babcock, Publications of the Astronomical Society of the Pacific **65**, 229 (1953).

[11] M. J. Booth, Light Sci Appl **3**, e165 (2014).

[12] M. Tanter, J.-F. Aubry, J. Gerber, J.-L. Thomas, and M. Fink, J. Acoust. Soc. Am. **110**, 37 (2001).

[13] A. Derode, A. Tourin, J. de Rosny, M. Tanter, S. Yon, and M. Fink, Phys. Rev. Lett. **90**, 014301 (2003).

[14] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Phys. Rev. Lett. **104**, 100601 (2010).

[15] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, Nat. Phys. **13**, 497 (2017).

[16] D. Bouchet, S. Rotter, and A. P. Mosk, Nat. Phys. **17**, 564 (2021).

[17] N. Bender, A. Yamilov, A. Goetschy, H. Yılmaz, C. W. Hsu, and H. Cao, Nat. Phys. **18**, 309 (2022).

[18] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Nat. Rev. Phys. **2**, 141 (2020).

[19] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Sci. Adv. **6**, eaay7170 (2020).
<figure><img src="image_24.png" /><figcaption>FIG. S18. Decomposition of the focusing law $\phi$ 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 bar: 50 µm). **c.** Associated confocal gain (scale bar: 50 µm). **d.** Example of local RPSF obtained in the white rectangle displayed in panel **b**1 (scale bar: 2 µm).</figcaption></figure>

sample and reference arms; (ii) an extremely long acquisition time since the focusing process
[20] S. Yoon, H. Lee, J. H. Hong, Y.-S. Lim, and W. Choi, Nat. Commun. **11**, 5721 (2020).

[21] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Proc. Natl. Acad. Sci. U. S. A. **117**, 14645 (2020).

[22] S. Gigan and et al., J. Phys. Photon. **4**, 042501 (2022), chapter 10.

[23] S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, Nat. Commun. **8**, 2157 (2017).

[24] Y. Kwon, J. H. Hong, S. Kang, H. Lee, Y. Jo, K. H. Kim, S. Yoon, and W. Choi, Nat. Commun. **14**, 105 (2023).

[25] A. Badon, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, Phys. Rev. Lett. **114**, 023901 (2015).

[26] A. Badon, D. Li, G. Lerosey, A. Claude Boccara, M. Fink, and A. Aubry, Optica **3**, 1160 (2016).

[27] E. Beaurepaire, A. C. Boccara, M. Lebec, L. Blanchot, and H. Saint-Jalmes, Opt. Lett. **23**, 244 (1998).

[28] A. Dubois, L. Vabre, A.-C. Boccara, and E. Beaurepaire, Appl. Opt. **41**, 805 (2002).

[29] V. Barolle, J. Scholler, P. Mecê, J.-M. Chassot, K. Groux, M. Fink, A. C. Boccara, and A. Aubry, Opt. Express **29**, 22044 (2021).

[30] J. Mertz, *Introduction to Optical Microscopy*, second edition ed. (Cambridge University Press, Cambridge, United Kingdom; New York, NY, 2019).

[31] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, Sci. Adv. **2**, e1600370 (2016).

[32] S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J.-S. Lee, Y.-S. Lim, Q.-H. Park, and W. Choi, Nat. Photonics **9**, 253 (2015).

[33] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry, Phys. Rev. X **10**, 021048 (2020).

[34] C. Brütt, A. Aubry, B. Gérardin, A. Derode, and C. Prada, Phys. Rev. E **106**, 025001 (2022).

[35] F. Roddier, ed., *Adaptive Optics in Astronomy* (Cambridge University Press, Cambridge, 1999).

[36] B. Judkewitz, R. Horstmeyer, I. M. Vellekoop, I. N. Papadopoulos, and C. Yang, Nat. Phys. **11**, 684 (2015).

[37] G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, Optica **4**, 886 (2017).
would have to be repeated 12 times on each of the 10⁸ points of the imaged volume. The
performance of matrix imaging is therefore impossible to reach with conventional adaptive
optics tools.

B. Comparison between iterative time reversal and iterative phase reversal

<figure><img src="image_25.png" /><figcaption>FIG. S19. 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. Scale bars: 50 µm.</figcaption></figure>

Even though the ITR process does not provide a satisfying estimator at the first iteration
of the matrix imaging process (see Supplementary Section S3 A and Supplementary Fig. S2),
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 Supplementary Figs. S19a\_2 and a\_3, we compare the confocal images provided
[38] W. Lambert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, IEEE Trans. Med. Imag. **41**, 3921 (2022).

[39] R. Mallart and M. Fink, J. Acoust. Soc. Am. **96**, 3721 (1994).

[40] S. Silverstein, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **48**, 1023 (2001).

[41] T. Wu, P. Berto, and M. Guillon, Appl. Phys. Lett. **118**, 251102 (2021).

[42] T. Wu, Y. Zhang, B. Blochet, P. Arjmand, P. Berto, and M. Guillon, arXiv:2304.01759v1 [cs.II], 10.48550/arXiv:2304.01759 (2023).

[43] K. Grieve, D. Ghoubay, C. Georgeon, G. Latour, A. Nahas, K. Plamann, C. Crotti, R. Bocheux, M. Borderie, T.-M. Nguyen, F. Andreiuolo, M.-C. Schanne-Klein, and V. Borderie, Sci. Rep. **7**, 13584 (2017).

[44] L. Zhu, J. B. de Monvel, P. Berto, S. Brasselet, S. Gigan, and M. Guillon, Optica **7**, 338 (2020).

[45] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, arXiv:2303.07483v1 [cs.II], 10.48550/arXiv:2303.07483 (2023).

[46] S. Kang, Y. Kwon, H. Lee, S. Kim, J. H. Hong, S. Yoon, and W. Choi, arXiv:2302.09503v1 [cs.II], 10.48550/arXiv:2302.09503 (2023).

[47] M. Jang, H. Ruan, I. M. Vellekoop, B. Judkewitz, E. Chung, and C. Yang, Biomed. Opt. Express **6**, 72 (2014).

[48] J. Scholler, K. Groux, O. Goureau, J.-A. Sahel, M. Fink, S. Reichman, C. Boccara, and K. Grieve, Light Sci. Appl. **9**, 140 (2020).

[49] B.-F. Osmanski, G. Montaldo, M. Tanter, and M. Fink, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **59**, 1575 (2012).

[50] M. Campillo and A. Paul, Science **299**, 547 (2003).

[51] A. Labeyrie, Astron. Astrophys. **6**, 85 (1970).

[52] S. Rotter and S. Gigan, Rev. Mod. Phys. **89**, 015005 (2017).

[53] G. Foschini and M. Gans, Wireless Personal Communications **6**, 311 (1998).

[54] A. F. Fercher, J. Biomed. Opt. **1**, 157 (1996).

[55] R. Bocheux, P. Pernot, V. Borderie, K. Plamann, and K. Irsch, PLOS ONE **14**, e0221707 (2019).

[56] U. Tricoli and R. Carminati, J. Opt. Soc. Am. A **36**, C122 (2019).

[57] A. Goicoechea, C. Brutt, F. Bureau, A. Le Ber, C. Prada, and A. Aubry, to be submitted (2023).
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 Supplementary Fig. S2 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.

To understand why, let us first express the confocal image as a function of the reflection matrix coefficients expressed in the plane wave basis:

$$I(\Delta\rho = 0, \rho_{\text{out}}, z) = \sum_{\mathbf{u}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) \exp \left[ i \frac{2\pi}{\lambda f} (\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}}) \cdot \rho_{\text{out}} \right] \quad (\text{S69})$$

Under an isoplanatic assumption, the coefficients of the reflection matrix in the plane wave basis can be expressed as follows²¹:

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

with $\tilde{\gamma}(\mathbf{u}, z) = \int d\rho\gamma(\rho, z) \exp(-i2\pi\mathbf{u}\cdot\rho/(\lambda f))$, the spatial frequency spectrum of the sample reflectivity at depth z. In absence of aberrations, each antidiagonal $(\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}} = \text{constant})$ encodes one spatial frequency of the sample reflectivity: $R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) = \tilde{\gamma}(\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}}, z)$ in the pupil support. $I(\Delta\rho = 0, \rho_{\text{out}}, z)$ is then a satisfying estimator of $\gamma(\rho, z)$ with a transverse resolution only limited by diffraction. In presence of aberrations, phase fluctuations of $\mathcal{T}_{\text{in/out}}$ implies phase distortions between each spatial frequency component of the object while amplitude variations of $\mathcal{T}_{\text{in/out}}$ implies the attenuation of some spatial frequency components of the object.

The compensation of aberrations consists in applying the phase conjugate of the aberration transmittance estimators, $\hat{T}_{\text{out}}$ and $\hat{T}_{\text{in}}$, at input and output of the reflection matrix:

$$R''(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) = \hat{T}_{\text{out}}^*(\mathbf{u}_{\text{out}}, z) R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) \hat{T}_{\text{in}}^*(\mathbf{u}_{\text{in}}, z) \quad (\text{S71})$$

Assuming ITR provides correct estimators of the aberration transmittance ($\hat{T}_{\text{in/out}} \equiv \mathcal{T}_{\text{in/out}}$) and injecting Eq. S70 into the last equation leads to the following expression for the reflection matrix corrected by ITR processing:

$$R''_{\text{ITR}}(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) = |\mathcal{T}_{\text{out}}(\mathbf{u}_{\text{out}}, z)|^2 \tilde{\gamma}(\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}}, z) |\mathcal{T}_{\text{in}}(\mathbf{u}_{\text{in}}, z)|^2 \quad (\text{S72})$$

Because the amplitude of $\mathcal{T}_{\text{in/out}}$ vanishes for high spatial frequency $|\mathbf{u}_{\text{in/out}}|$ (Supplementary Fig. S2b), ITR tends to filter the high spatial frequency component of the object.
[58] J.-L. Robert and M. Fink, J. Acoust. Soc. Am. **123**, 866 (2008).

[59] C. Prada and M. Fink, Wave Motion **20**, 151 (1994).

[60] C. Prada, S. Manneville, D. Spoliansky, and M. Fink, J. Acoust. Soc. Am. **99**, 2067 (1996).
On the contrary, IPR converges towards a normalized version of the transmittances $T_{in/out}$, such that: $\hat{T}_{in/out} = \exp(j\arg\{T_{in/out}\})$. Using this last expression of $\hat{T}_{in/out}$ and injecting Eq. S70 into S71 leads to the following expression for the reflection matrix corrected by IPR processing:

$$R''_{IPR}(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, z) = |T_{\text{out}}(\mathbf{u}_{\text{out}}, z)|\tilde{\gamma}(\mathbf{u}_{\text{in}} + \mathbf{u}_{\text{out}}, z)|T_{\text{in}}(\mathbf{u}_{\text{in}}, z)| \quad (S73)$$

IPR leaves the spatial frequency spectrum of the reflection matrix unchanged and does not filter the high spatial frequency components of the reflection matrix. Of course, a better compensation of aberration would consist in an inverse filter where the amplitude decrease of $T_{in/out}$ could be compensated on top of phase distortions. Nevertheless, this operation is extremely sensitive to noise and suppose a perfect match between the aberration transmittance $T_{in/out}$ and their estimators $\hat{T}_{in/out}$. This is wrong especially for the phase of $T_{in/out}$ at high spatial frequencies where the amplitude of $\hat{T}_{in/out}$ vanishes. Therefore, IPR is an adequate compromise between the matched filter operated by ITR that affect the high spatial frequencies of the object and an inverse filter that is extremely sensitive to noise.

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. Supplementary Fig. S19a_1 shows the ratio between their confocal gains at depth $z=200~\mu 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~\mu m$), the ITR process leads to an image with strong vignetting effects (Supplementary Fig. S19b_3), which is a manifestation of a lack of correlations of the estimator $\hat{T}$ between adjacent windows. On the contrary, the IPR leads to an estimator whose scattered component $\hat{S}$ exhibits the expected memory effect for this size of spatial windows (Supplementary Fig. S7d). Supplementary Figure S19 thus shows the superiority of IPR in terms of spatial resolution for the $T$-matrix estimation.

The superiority of IPR compared with ITR is also confirmed by the resolution target
experiment presented in Supplementary Section S3 H. Supplementary Figure S20 compares the images obtained via IPR and ITR using the multi-scale process described above. The full-field images show a slightly better contrast with IPR (Supplementary Fig. S20b) than ITR (Supplementary Fig. S20c). The zoom on the smallest patterns of the target also shows the benefit of IPR (Supplementary Fig. S20f) compared with ITR (Supplementary Fig. S20g) in terms of resolution. Indeed, the three strips are clearly resolved with IPR, while this pattern remains quite blurred with ITR. This superiority is also confirmed by the corresponding RPSFs displayed in Supplementary Figs. S20j and k. In the area of the smallest pattern surrounded by a white circle, the final RPSFs obtained with IPR display a weaker incoherent background than with ITR. This reference experiment thus confirms the overall superiority of IPR compared with ITR that we have already noticed in the cornea experiment.

### C. 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 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. Supplementary Figure S21 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 [Supplementary Figs. S21b$_2$ and c$_2$], the obtained images show significant differences [Supplementary Figs. S21b$_1$ and c$_1$]. On the one hand, the
<figure><img src="image_26.png" /><figcaption>FIG. S20. Comparison of IPR with other state-of-the-art methods for the resolution target experiment. a-d. Confocal images (scale bar: 50 µm): Original image (a), corrected image using the IPR multi-scale process (b), corrected image using the ITR multi-scale process (c) and corrected image using the local CLASS algorithm (d). e-h. Corresponding blowups of the area contained inside the blue rectangle displayed in panels a-d. i-k. Corresponding maps of RPSFs (computed over spatial windows of size L = 2 µm).</figcaption></figure>

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 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 Supplementary Fig. S22 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 Supplementary Fig. S22b<sub>2</sub>. It shows an important background, a manifestation of an imperfect compensation of forward multiple
<figure><img src="image_27.png" /><figcaption>FIG. S21. Comparison between CLASS and RMI upon convergence. **a.** Original confocal image (scale bar: 50 µm) (**a**₁) and mean RPSF (scale bar: 2 µm) (**a**₂) at z=200 µm. **b-c.** Corresponding CLASS (**b**₁) and IPR (**c**₁) images with their mean RPSFs (**b**₂ and **c**₂, respectively). The size of spatial windows is *L* = 13 µm.</figcaption></figure>

scattering. On the contrary, the RPSF obtained at the same location via IPR is close to be ideal (Supplementary Fig. S22c₂). When looking into details to the CLASS image (Supplementary Fig. S22b₁), strong vignetting effects can be observed while the IPR image exhibits a continuous reflectivity (Supplementary Fig. S22c₁).

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). Supplementary Figure S23 shows this correlation map for the mid-point **r**<sub>p</sub> of the area displayed with a white square in Supplementary Fig. S22. While the **S**-matrix derived by IPR preserves a short-range correlation between neighbouring windows (see Supplementary Fig. S23c), the CLASS algorithm leads to a fully spatially incoherent estimator **S**̂ (see Supplementary Fig. S23a). 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 Supplementary Fig. S23 accounts for the vignetting effects observed for the CLASS image in Supplementary Fig. S22b₁.
<figure><img src="image_28.png" /><figcaption>FIG. S22. Comparison between CLASS and RMI for small isoplanatic patches. **a.** Original confocal image (scale bar: 50 µm) (**a**₁) and mean RPSF (scale bar: 2 µm) (**a**₂) at *z* =200 µm. **b-c.** Corresponding CLASS (**b**₁) and IPR (**c**₁) images with their mean RPSFs (**b**₂ and **c**₂, respectively). The size of spatial windows is *L* = 6 µm.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIG. S23. Memory effect exhibited the forward multiple scattering estimator: Comparison between CLASS and RMI. Spatial correlation function $C_S(\mathbf{r}_p, \mathbf{r}'_p)$ of the scattering component $\hat{\mathbf{S}}$ for a given point $\mathbf{r}'_p$. **a.** CLASS. **b.** RMI. The depth is $z = 200$ µm and the size of spatial windows is $L = 6$ µm. Scale bar: 50 µm.</figcaption></figure>

This vignetting effect is also highlighted by the resolution target experiment in Supplement
mentary Fig. S20. The corresponding images obtained with the multi-scale IPR process and
the local CLASS algorithm are displayed in Supplementary Figs. S20b and d, respectively.
For each method, the size of spatial windows over which aberration phase laws are estimated
is *L* = 10 µm. As in the cornea experiment (Supplementary Fig. S23), the CLASS image
lacks of continuity between adjacent spatial windows (Supplementary Fig. S20d). The map of
RPSFs also shows the imperfect convergence of the CLASS process (Supplementary Fig. S20l).
Unlike the RPSFs obtained at the end of the IPR process (Supplementary Figs. S20j), the
CLASS RPSFs are far from displaying a diffraction-limited feature. They also show a strong
variability across the field-of-view, which is another manifestation of the vignetting effect
observed in the full-field CLASS image (Supplementary Figs. S20d). The resolution target
experiment thus confirms the superiority of a multi-scale IPR process with respect to the
local CLASS algorithm.
S6. GLOSSARY

<table><caption>TABLE S1. Glossary of the operators used in this study</caption><thead><tr><th>operator</th><th>definition</th></tr></thead><tbody><tr><td><b>R</b> = [R(ρ<sub>in</sub>, ρ<sub>out</sub>, z)]</td><td>reflection matrix in the focused basis</td></tr><tr><td><b>R</b><sub>in/out</sub> = [R<sub>in</sub>(Δρ<sub>in/out</sub>, ρ<sub>out/in</sub>, z)]</td><td>de-scanned matrix at input/output</td></tr><tr><td><b>D</b><sub>in/out</sub> = [D(<b>u</b><sub>in/out</sub>, ρ<sub>out/in</sub>), z]</td><td>input/output distortion matrix</td></tr><tr><td><b>H</b><sub>in/out</sub> = [H<sub>in</sub>(ρ<sub>s</sub>, ρ<sub>in/out</sub>, z)]</td><td>input/output focusing matrix</td></tr><tr><td><b>T</b><sub>0</sub> = [T<sub>0</sub>(<b>u</b>, ρ)]</td><td>Fourier transform operator</td></tr><tr><td><b>C</b><sub>in/out</sub> = [C<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>u'</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output pupil correlation matrix of <b>D</b><sub>in</sub></td></tr><tr><td><b>C</b><sub>H</sub> = [C<sub>H</sub>(<b>u</b>, <b>u'</b>)]</td><td>correlation matrix associated with<br>a virtual reflector of reflectivity |H(ρ<sub>s</sub>)|^2</td></tr><tr><td><b>T</b><sub>in/out</sub> = [T<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output transmission matrix</td></tr><tr><td><b>T</b><sub>in/out</sub> = [<b>T</b><sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output transmittance matrix</td></tr><tr><td><b>T̂</b><sub>in/out</sub> = [<b>T</b><sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>estimator of the input/output transmittance matrix</td></tr><tr><td>ϕ<sub>in/out</sub> = [ϕ<sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>phase of the input/output transmission matrix estimator</td></tr><tr><td><b>A</b><sub>in/out</sub> = [<b>A</b><sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output aberration matrix</td></tr><tr><td><b>S</b><sub>in/out</sub> = [<b>S</b><sub>in/out</sub>(<b>u</b><sub>in/out</sub>, <b>r</b><sub>p</sub>)]</td><td>input/output forward multiple scattering matrix</td></tr></tbody></table>
References.

[1] B. Gérardin, J. Laurent, A. Derode, C. Prada, and A. Aubry, Full transmission and reflection of waves propagating through a maze of disorder, Phys. Rev. Lett. **113**, 173901 (2014).

[2] M. Horodynski, M. Kuhmayer, C. Ferise, S. Rotter, and M. Davy, Anti-reflection structure for perfect transmission through complex media, Nature **607**, 281 (2022).

[3] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Direct observation of anderson localization of matter waves in a controlled disorder, Nature **453**, 891 (2008).

[4] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Localization of ultrasound in a three-dimensional elastic network, Nat. Phys. **4**, 945 (2008).

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

[6] J. Bertolotti and O. Katz, Imaging in complex media, Nat. Phys. **18**, 1008 (2022).

[7] C. Dunsby and P. French, Techniques for depth-resolved imaging through turbid media including coherence-gated imaging, J. Phys. D: Appl. Phys. **36**, R207 (2003).

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

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

[10] H. W. Babcock, The possibility of compensating astronomical seeing, Publications of the Astronomical Society of the Pacific **65**, 229 (1953).

[11] M. J. Booth, Adaptive optical microscopy: The ongoing quest for a perfect image, Light Sci Appl **3**, e165 (2014).

[12] M. Tanter, J.-F. Aubry, J. Gerber, J.-L. Thomas, and M. Fink, Optimal focusing by spatio-temporal inverse filter. i. basic principles, J. Acoust. Soc. Am. **110**, 37 (2001).

[13] A. Derode, A. Tourin, J. de Rosny, M. Tanter, S. Yon, and M. Fink, Taking advantage of multiple scattering to communicate with time-reversal antennas, Phys. Rev. Lett. **90**, 014301 (2003).

[14] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media, Phys. Rev. Lett. **104**, 100601 (2010).
[15] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, Correlation-enhanced control of wave focusing in disordered media, Nat. Phys. **13**, 497 (2017).

[16] D. Bouchet, S. Rotter, and A. P. Mosk, Maximum information states for coherent scattering measurements, Nat. Phys. **17**, 564 (2021).

[17] N. Bender, A. Yamilov, A. Goetschy, H. Yılmaz, C. W. Hsu, and H. Cao, Depth-targeted energy delivery deep inside scattering media, Nat. Phys. **18**, 309 (2022).

[18] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Deep optical imaging within complex scattering media, Nat. Rev. Phys. **2**, 141 (2020).

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

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

[21] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media, Proc. Natl. Acad. Sci. U. S. A. **117**, 14645 (2020).

[22] S. Gigan and et al., Roadmap on wavefront shaping and deep imaging in complex media, J. Phys. Photon. **4**, 042501 (2022), chapter 10.

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

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

[25] A. Badon, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, Retrieving Time-Dependent Green's Functions in Optics with Low-Coherence Interferometry, Phys. Rev. Lett. **114**, 023901 (2015).

[26] A. Badon, D. Li, G. Lerosey, A. Claude Boccara, M. Fink, and A. Aubry, Spatio-temporal imaging of light transport in highly scattering media under white light illumination, Optica **3**, 1160 (2016).
[27] E. Beaurepaire, A. C. Boccara, M. Lebec, L. Blanchot, and H. Saint-Jalmes, Full-field optical coherence microscopy, Opt. Lett. **23**, 244 (1998).

[28] A. Dubois, L. Vabre, A.-C. Boccara, and E. Beaurepaire, High-resolution full-field optical coherence tomography with a Linnik microscope, Appl. Opt. **41**, 805 (2002).

[29] V. Barolle, J. Scholler, P. Mecê, J.-M. Chassot, K. Groux, M. Fink, A. C. Boccara, and A. Aubry, Manifestation of aberrations in full-field optical coherence tomography, Opt. Express **29**, 22044 (2021).

[30] J. Mertz, *Introduction to Optical Microscopy*, second edition ed. (Cambridge University Press, Cambridge, United Kingdom; New York, NY, 2019).

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

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

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

[34] M. P. V. Albada and A. Lagendijk, Observation of weak localization of light in a random medium, Phys. Rev. Lett. **55**, 2692 (1985).

[35] P.-E. Wolf and G. Maret, Weak localization and coherent backscattering of photons in disordered media, Phys. Rev. Lett. **55**, 2696 (1985).

[36] A. Goicoechea, C. Brutt, F. Bureau, A. Le Ber, C. Prada, and A. Aubry, Reflection measurement of the scattering mean free path at the onset of multiple scattering, arXiv:2401.02246 10.48550/arXiv.2401.02246 (2024).

[37] F. Roddier, ed., *Adaptive Optics in Astronomy* (Cambridge University Press, Cambridge, 1999).

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

[39] G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, Generalized optical memory effect, Optica **4**, 886 (2017).

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

[41] R. Mallart and M. Fink, Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion, J. Acoust. Soc. Am. **96**, 3721 (1994).

[42] T. Wu, Y. Zhang, B. Blochet, P. Arjmand, P. Berto, and M. Guillon, Single-shot digital optical fluorescence phase conjugation through forward multiple-scattering samples, Sci. Adv. **10**, ad1120 (2023).

[43] K. Grieve, D. Ghoubay, C. Georgeon, G. Latour, A. Nahas, K. Plamann, C. Crotti, R. Bocheux, M. Borderie, T.-M. Nguyen, F. Andreiuolo, M.-C. Schanne-Klein, and V. Borderie, Stromal striae: a new insight into corneal physiology and mechanics, Sci. Rep. **7**, 13584 (2017).

[44] L. Zhu, J. B. de Monvel, P. Berto, S. Brasselet, S. Gigan, and M. Guillon, Chromato-axial memory effect through a forward-scattering slab, Optica **7**, 338 (2020).

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

[46] S. Kang, Y. Kwon, H. Lee, S. Kim, J. H. Hong, S. Yoon, and W. Choi, Tracing multiple scattering trajectories for deep optical imaging in scattering media, Nat. Commun. **14**, 6871 (2023).

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

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

[49] B.-F. Osmanski, G. Montaldo, M. Tanter, and M. Fink, Aberration correction by time reversal of moving speckle noise, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **59**, 1575 (2012).

[50] M. Campillo and A. Paul, Long-range correlations in the diffuse seismic coda, Science **299**, 547 (2003).

[51] A. Labeyrie, Attainment of diffraction limited resolution in large telescopes by fourier analysing speckle patterns in star images, Astron. Astrophys. **6**, 85 (1970).

[52] S. Rotter and S. Gigan, Light fields in complex media: Mesoscopic scattering meets wave control, Rev. Mod. Phys. **89**, 015005 (2017).
[53] G. Foschini and M. Gans, On limits of wireless communications in a fading environment when using multiple antennas, Wireless Personal Communications **6**, 311 (1998).

[54] U. Najar, V. Barolle, P. Balondrade, M. Fink, C. Boccara, and A. Aubry, Optical distortion matrices deep inside a human cornea [data]. Zenodo (2023).

[55] A. F. Fercher, Optical coherence tomography., J. Biomed. Opt. **1**, 157 (1996).

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

[57] C. Prada and M. Fink, Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media, Wave Motion **20**, 151 (1994).

[58] C. Prada, S. Manneville, D. Spoliansky, and M. Fink, Decomposition of the time reversal operator: Detection and selective focusing on two scatterers, J. Acoust. Soc. Am. **99**, 2067 (1996), publisher: Acoustical Society of America.

[59] S. Silverstein, Ultrasound scattering model: 2-d cross-correlation and focusing criteria-theory, simulations, and experiments, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **48**, 1023 (2001).

[60] T. Wu, P. Berto, and M. Guillon, Appl. Phys. Lett. **118**, 251102 (2021).

[61] E. Akkermans, P. E. Wolf, R. Maynard, and G. Maret, Theoretical study of the coherent backscattering of light by disordered media, J. Phys. France **49** (1988).

[62] A. Tourin, A. Derode, P. Roux, B. A. van Tiggelen, and M. Fink, Time-dependent coherent backscattering of acoustic waves, Phys. Rev. Lett. **79**, 3637 (1997).

[63] E. Larose, L. Margerin, B. A. van Tiggelen, and M. Campillo, Weak localization of seismic waves, Phys. Rev. Lett. **93**, 048501 (2004).

[64] A. Aubry and A. Derode, Ultrasonic imaging of highly scattering media from local measurements of the diffusion constant: Separation of coherent and incoherent intensities, Phys. Rev. E **75**, 026602 (2007).

[65] L. A. Cobus, G. Maret, and A. Aubry, Crossover from renormalized to conventional diffusion near the 3D anderson localization transition for light, Phys. Rev. B **106**, 014208 (2022).

[66] W. Lambert, J. Robin, L. A. Cobus, M. Fink, and A. Aubry, Ultrasound matrix imaging – Part I: The focused reflection matrix, the F-factor and the role of multiple scattering, IEEE Trans. Med. Imag. **41**, 3907 (2022).
[67] R. Bocheux, P. Pernot, V. Borderie, K. Plamann, and K. Irsch, Quantitative measures of corneal transparency, derived from objective analysis of depth-resolved corneal images, demonstrated with full-field optical coherence tomographic microscopy, PLOS ONE 14, e0221707 (2019).

[68] U. Tricoli and R. Carminati, Modeling of full-field optical coherence tomography in scattering media, J. Opt. Soc. Am. A 36, C122 (2019).

[69] C. Brütt, A. Aubry, B. Gérardin, A. Derode, and C. Prada, Weight of single and recurrent scattering in the reflection matrix of complex media, Phys. Rev. E 106, 025001 (2022).