Smart optical coherence tomography for ultradeep imaging through highly scattering media
Abstract
Multiple scattering of waves in disordered media is a nightmare whether it be for detection or imaging purposes. The best approach so far to get rid of multiple scattering is optical coherence tomography. It basically combines confocal microscopy and coherence timegating to discriminate ballistic photons from a predominant multiple scattering background. Nevertheless, the imaging depth range remains limited to 1 mm at best in human soft tissues. Here we propose a matrix approach of optical imaging to push back this fundamental limit. By combining a matrix discrimination of ballistic waves and iterative timereversal, we show both theoretically and experimentally an extension of the imagingdepth limit by at least a factor two compared to optical coherence tomography. In particular, the reported experiment demonstrates imaging through a strongly scattering layer from which only one reflected photon over 1000 billion is ballistic. This approach opens a new route towards ultradeep tissue imaging.
Introduction
The propagation of light in inhomogeneous media is a fundamental problem with important applications, ranging from astronomical observations through a turbulent atmosphere to deep tissue imaging in microscopy or light detection through a dense cloud in LIDAR technology. Conventional focusing and imaging techniques based on the Born approximation generally fail in strongly scattering media due to the multiple scattering (MS) events undergone by the incident wavefront. Recent advances in light manipulation techniques have allowed great progresses in optical focusing through complex media Mosk_review (). Following pioneering works in ultrasound derode0 (); derode_2003 (), Vellekoop and Mosk Vellekoop () showed how light can be focused spatially through a strongly scattering medium by actively shaping the incident wavefront with a spatial light modulator (SLM). Subsequently, a matrix approach of light propagation through complex media was developed popoff (). It relies on the measurement of the Green’s functions between each pixel of a SLM and of a chargecoupled device (CCD) camera across a scattering medium. The experimental access to this transmission matrix allows to take advantage of MS for optimal light focusing popoff (); choi (); salma () and communication popoff2 (); Cizmar () across a diffusive layer or a multimode fiber.
MS is a much more difficult challenge with regards to imaging. On the one hand, imaging techniques like diffuse optical tomography Durduran (), acoustooptic Resink () or photoacoustic Xu () imaging take advantage of the diffuse light to image scattering media in depth but their resolution power is limited. More recently, the memoryeffect exhibited by the MS speckle freund (); feng () has been taken advantage of to image objects through strongly scattering layers with a diffractionlimited resolution bertolotti (); katz_2012 (); katz_2014 (). However, it only applies to thin opaque layers as the fieldofview is inversely proportional to the scattering medium thickness. On the other hand, conventional reflection imaging methods provide an optical diffractionlimited resolution but usually rely on a single scattering (SS) assumption. The imaging limit of conventional microscopy can be derived from the scattering mean free path . It describes the average distance that a photon travels between two consecutive scattering events. In turbid media, MS starts to predominate beyond a few . To cope with the fundamental issue of MS, several approaches have been proposed in order to enhance the SS contribution drowned into a predominant MS background. The first option is to spatially discriminate SS and MS as performed in confocal microscopy Pawley (); Ntziachristos () or twophoton microscopy Theer (). The second option consists in separating SS from MS photons by means of time gating Fujimoto2 (); Alfano (); Alfano2 (). Probably, the most widely employed coherent timegated technique is optical coherence tomography (OCT) fujimoto (); fujimoto3 (); Fercher (), which is the analogous to ultrasound imaging. It combines scanning confocal microscopy with coherent heterodyne detection oct0 (). OCT has drastically extended the imagingdepth limit compared to conventional microscopy. Nevertheless, its ability of imaging soft tissues remains typically restricted to a depth of 1 mm oct (); Dunsby ().
Inspired by previous works in ultrasound imaging through strongly scattering media aubry (); aubry2 (); shahjahan (), we propose a matrix approach of optical imaging to push back the fundamental limit of MS. Experimentally, this approach, referred to as smartOCT, relies on the measurement of a timegated reflection matrix from the scattering medium. Unlike previous works choi2 (); Kang (), the reflection matrix is here directly investigated in the focal plane on a pointtopoint basis robert (); robert2 (). An inputoutput analysis of the reflection matrix allows to get rid of most of the MS contribution. Iterative timereversal prada (); prada2 (); popoff3 () is then applied to overcome the residual MS contribution as well as the aberration effects induced by the turbid medium itself. A proofofconcept experiment demonstrates imaging through a strongly scattering paper layer from which only one reflected photon over 1000 billion is associated to a SS event from the object hidden behind it. In a second experiment, our approach is successfully applied to optical imaging through a thick layer of biological tissues. Compared to OCT and related methods Kang (), we show both theoretically and experimentally an extension of the imagingdepth limit by a factor two. This means a multiplication by a factor two of the sensitivity in dB compared to existing OCT systems.
Results
Measuring the timegated reflection matrix
The smartOCT approach is based on the measurement of a timegated reflection matrix from the scattering sample. Until now, optical transmission/reflection matrices have always been measured in the space (plane wave basis) popoff (); choi (); choi2 (); Kang (). Here, inspired by previous studies in acoustics robert (); robert2 (), the reflection matrix is directly investigated in the real space (pointtopoint basis). The experimental set up and procedure are described in Fig. 1. A set of reflection coefficients are measured between each point of the focal plane identified by the vectors at the input and at the output. These coefficients form the reflection matrix . In the experiments described below, the reflection matrix is measured over a fieldofview of m mapped with input wavefronts.
Reflection matrix in the single scattering regime
Figure 2A displays a reference reflection matrix measured for a ZnO bead of diameter m deposited on a microscope slide. contains two main contributions:

The specular echo from the glass slide that emerges throughout the diagonal of

The strong bead echo that arises for positions , with the position of the bead in the focal plane. Because the bead diameter is larger than the width of the point spread function, the bead also emerges along offdiagonal elements for which .
The single scattering contribution thus only emerges along the diagonal and closeddiagonal elements of . This is accounted for by the fact that a singlyscattered wavefield can only come from points illuminated by the incident focal spot. A timegated confocal image can be deduced from by only considering its diagonal elements, such that
(1) 
The corresponding image displayed in Fig. 2B is equivalent to a enface OCT image. Not surprisingly, it shows a clear image of the target on the microscope glass slide.
Reflection matrix in the deep multiple scattering regime
In a second experiment, a stack of two paper sheets is placed between the MO and the target bead [see Fig. 1]. The thickness of each sheet is of 82 m, hence an overall thickness of m for the scattering layer. The distance between the frontsurface of the scattering layer and the target is mm. The scattering and transport mean free paths in the paper sheet have been measured and estimated to be m and m, respectively [see Supplementary Section I]. This yields an optical thickness . The ballistic wave has to go through 24.5 back and forth, thus undergoing an attenuation of in intensity. The singletomultiple scattering ratio (SMR) of the reflected wavefield is estimated to be in the MO backfocal plane [see Supplementary Section II]. For an incident plane wave, it means that only 1 scattered photon over 1000 billion is associated to a SS event from the target. As shown in Supplementary Fig. S3, the target is far to be detectable and imaged in this experimental configuration, whether it be by conventional microscopy (SMR ), confocal microscopy (SMR ) or by OCT (SMR ). This experimental situation is thus particularly extreme, even almost desperate, for a successful imaging of the target.
Figure 2C displays the reflection matrix measured in presence of the scattering layer. Contrary to the SS contribution that emerges along diagonal and closeddiagonal elements of [Fig. 2A], MS randomizes the directions of light propagation and gives rise to a random reflection matrix aubry (). Nevertheless, one can try to image the target by considering the diagonal of [equation (1)]. The corresponding enface OCT image is shown in Fig. 2D. As theoretically expected [see Supplementary Fig. S3], MS still predominates despite confocal filtering and coherence timegating. An image of speckle is thus obtained without any enhancement of the intensity at the expected target location [see the comparison with the reference image in Fig. 2B].
Matrix approach dedicated to target detection in the deep multiple scattering regime
An alternative route is now proposed to image and detect the target behind the scattering layer. The smartOCT approach first consists in filtering the multiplescattering contribution in the measured reflection matrix . To that aim, the matrix is projected on a characteristic SS matrix , whose elements are given by
(2) 
is a tunable parametric length that accounts for the fact that the ballistic signal does not only emerge along the diagonal of the matrix but also along offdiagonal elements [Fig. 2A]. is governed by two factors:

The coherence length of the ballistic wavefield in the focal plane: In addition to ballistic attenuation and MS, the scattering layer also induces aberrations that degrades the focusing quality of the ballistic wavefront and enlarge the point spread function of the imaging system.

The size of the target: As shown by Fig. 2A, the target signal does not only emerge along the diagonal elements of in absence of the scattering layer. This is accounted for by the size of the target which is larger than the resolution cell.
Mathematically, the projection of can be expressed as an Hadamard product with ,
(3) 
which, in term of matrix coefficients, can be written as
(4) 
This mathematical operation thus consists in keeping the diagonal and closeddiagonal coefficients of where the SS contribution arises and filtering the offdiagonal elements of mainly associated with the MS contribution. It can be seen as a digital confocal operation with a virtual pinhole mask of size psaltis (). In the present experiment, a SS matrix is deduced from by considering m. The result is displayed in Fig. 2E. contains the SS contribution as wanted plus a residual MS contribution [see the comparison with the reference matrix in Fig. 2A]. This term persists because MS signals also arise along and close to the diagonal of . Compared to a single/multiple scattering separation performed in the farfield aubry2 (); Kang (), a single scattering projection in a pointtopoint basis [Eq.3] is much more flexible since the tunable parameter can be adapted as a function of the aberration level or the expected target size.
Once this SS matrix is obtained, one can apply the DORT method (French acronym for Decomposition of the Time Reversal Operator). Initially developed for ultrasound prada (); prada2 (), the DORT method takes advantage of the reflection matrix to focus iteratively by time reversal processing on each scatterer of a multitarget medium popoff3 (). Mathematically, the timereversal invariants can be deduced from the eigenvalue decomposition of the time reversal operator or, equivalently, from the singular value decomposition of (the superscript stands for transpose conjugate). A onetoone association between each eigenstate of and each scatterer does exist. On the one hand, the eigenvectors of allow selective focusing and imaging of each scatterer. On the other hand, the associated eigenvalue directly yields the scatterer reflectivity. Nevertheless, this onetoone association is only valid under a SS approximation. Hence the DORT method cannot be applied to the raw matrix since it contains an extremely predominant MS contribution. The trick here is to take advantage of the SS matrix .
A singular value decomposition (SVD) of is performed. It consists in writing . is a diagonal matrix containing the real positive singular values in a decreasing order . The square of the singular values, , correspond to the eigenvalues of . and are unitary matrices whose columns correspond to the input and output singular vectors and , respectively. Fig. 2F displays the histogram of the eigenvalues normalized by their average. It is compared to the distribution that would be obtained in a fully multiple scattering regime [see Supplementary Section III]. The histogram of in Fig. 2F follows this distribution except for the largest eigenvalue . The latter one is actually beyond the superior bound of the MS continuum of eigenvalues. This means that the first eigenspace is associated to the target aubry2 (); shahjahan (). The combination of the first input and output singular vectors, , forms the smartOCT image displayed in Fig. 2G. The image of the target is nicely recovered. The comparison with the enface OCT image displayed in Fig. 2D unambiguously demonstrates the benefit of smartOCT in detecting a target in the deep MS regime ( ). Note that the target image does not match exactly with the reference image [Fig. 2B]. This difference can be accounted for by residual aberration effects induced by the scattering layer itself.
The theoretical study developed in the Supplementary Section II confirms this experimental result. In the conditions of the reported experiment, the imaging thicknesslimit (SMR) is found to be of 1.5 for conventional microscopy, 3.5 for confocal microscopy and 7 for OCT. This explains the failure of OCT in detecting the target [Fig. 2D]. On the contrary, the predicted imagingthickness limit is of 12 for the smartOCT approach. This accounts for the successful detection of the target in our experimental configuration [Fig. 2G].
Imaging in the deep multiple scattering regime
From the previous experiment, one could say that the smartOCT approach is only a single target detection method dedicated to the deep MS regime. To demonstrate that we can go beyond the detection and imaging of a single target, a configuration with three 5mdiameter ZnO beads deposited on a microscope slide is investigated. Figure 3 (A and B) display the timegated reflection matrix and the corresponding enface OCT image in absence of any scattering layer. Both figures will serve as reference in the following. A paper sheet is then placed between the MO and the focal plane. The corresponding reflection matrix is shown in Fig. 3C. Not surprisingly, it displays a random feature characteristic of a predominant MS contribution. The enface OCT image built from the diagonal elements of is displayed in Fig. 3D. Again the MS speckle prevents from a clear and unambiguous detection of the three targets [see the comparison with Fig. 3B]. The MS filter is applied to the raw matrix [equation (3), m], yielding a SS matrix displayed in Fig. 3E. Iterative timereversal is then performed. The three first eigenstates of are displayed in Fig. 3F. A comparison with the reference image in Fig. 3B highlights the onetoone association between each bead and each of these eigenstates. The combination of these eigenstates weighted by the corresponding eigenvalue finally provide the smartOCT image. The comparison with the enface OCT image [Fig. 3D] demonstrates the success of the smartOCT approach for imaging in the deep MS regime.
Imaging through thick biological tissues
Following this experimental proofofconcept, we now apply our approach to the imaging of an extended object through biological tissues. A positive USAF 1951 resolution target placed behind a 800 mthick layer of rat intestine tissues is imaged through an immersion objective (Olympus, 40, NA=0.8) [see Fig. 4A]. The reflection matrix is measured over a fieldofview of m [see the green square in Fig. 4A] with input wavefronts. The diagonal of yields the enface OCT image displayed in Fig. 4B. Due to the aberration effects and multiple scattering events induced by the biological tissues, the three bars of the USAF target cannot be recovered. A comparable result is obtained if the DORT method is directly applied to the raw matrix [See Supplementary Fig. S4]. To overcome these detrimental effects, the MS filter is applied to the raw matrix [equation (3), m], yielding a SS matrix . Iterative timereversal is then performed. In previous experiments, the object to image consisted in one or a few beads. This sparsity implied that only few eigenstates were needed to recover the image of the beads. In the present case, the USAF target is an extended object. It is thus associated with a large number of eigenstates, scaling as the number of resolution cells contained in the object robert3 (). To estimate the rank of the object, one can compute the standard deviation of the image, , as a function of the number of eigenstates considered for the imaging process. The result is displayed in Fig. 4C. A maximum standard deviation is found for eigenstates. The corresponding image is displayed in Fig. 4E. The three bars of the USAF target are nicely recovered and the comparison with the enface OCT image [Fig. 4B] is striking. This experimental result demonstrates the benefit of our approach for deeptissue imaging. To illustrate the importance of a correct determination of , we also show for comparison the images built from the 20 and 500 first eigenstates of [see Fig. 4(D and F), respectively]. On the one hand, considering too few eigenstates only provides a partial imaging of the field of view [Fig. 4D]. On the other hand, considering too many eigenstates blurs the image since the weakest eigenvalues are mainly associated with the multiple scattering background [Fig. 4F].
Discussion
These experimental results demonstrate the benefit of the smartOCT approach for optical imaging in strongly scattering media. In the Supplementary Section II, this superiority is confirmed by the theoretical investigation of the imagingdepth limit derived for different imaging techniques (conventional/confocal microscopy, OCT, smartOCT). In view of applications to deep tissue imaging, Fig. 5 shows the SMR evolution versus the optical depth in biological tissues. The experimental parameters chosen for this figure are those typically encountered in fullfield OCT Assayag () and are provided in Supplementary Tab.S1. The detection threshold is set at a SMR of 1. The imagingdepth limit expected in tissues is of 1 for conventional microscopy, 8 for confocal microscopy and 12 for OCT. The latter value is in agreement with the imagingdepth limit recently reached by Kang et al. Kang () with an optical technique similar to OCT. It actually combines coherence time gating with a spatial inputoutput correlation of waves from the farfield that allows a confocal discrimination of reflected photons. On the contrary, our approach goes beyond OCT as it also involves a subsequent iterative timereversal processing of the reflection matrix. It results in an additional gain in SMR that scales with , the number of input wavefronts (see Supplementary Section II). This leads to an imagingdepth limit of 22 in Fig. 5, hence multiplying by almost two the current OCT limit. Such an imaging improvement is drastic if we keep in mind that the ballistic contribution decreases by a factor in a reflection configuration. The smartOCT approach is thus particularly suited for ultradeep tissue imaging. Of course, a tradeoff will have to be made in practice between the imaging depth and the measurement time that also scales linearly with . Note also that the imaging depth can be limited by the dynamic range of the CCD camera. For instance, a dynamic range of dB would be required to reach the theoretical imaging depthlimit in the experimental conditions of Fig. 5 (see Supplementary Tab. S1).
A second point we would like to discuss is the resolution and sectioning capabilities of our approach. On the one hand, as the image is built from singlyscattered photons, the transverse resolution is only diffractionlimited and does not depend on the penetration depth. On the other hand, the coherent timegated detection scheme provides an axial resolution that is only governed by the coherence time of the light source: , with the refractive index of the medium. In the present experiment, the coherence time of the femtosecond laser is of 50 fs. It yields an axial resolution of 5 m in biological tissues (). By measuring a set of reflection matrices at successive depths, a 3D image of the sample can thus be obtained as in conventional OCT with the great advantage that the penetration depth is multiplied by a factor 2. To reach a better axial resolution, the measurement of the reflection matrix can be made under a simple white light illumination. A recent study has actually demonstrated the passive measurement of the timedependent pointtopoint reflection matrix from an incoherent illumination badon (). As already demonstrated in fullfield OCT Assayag (), the temporal incoherence of the whitelight source would provide an excellent axial resolution (1 m). Moreover, such a device would comply with the noninvasive, lowcost, speed and lowcomplexity specifications required for medical applications.
In summary, this study proposes a matrix approach of light propagation dedicated to optical detection and imaging through complex media. The socalled smartOCT approach combines a matrix discrimination of ballistic waves with iterative timereversal. A first proofofconcept experiment demonstrates the imaging of several microbeads in the deep multiple scattering regime, whereas existing imaging techniques such as OCT are shown to fail. A second experiment demonstrates the diffractionlimited imaging of an extended object (USAF target) through a thick layer of biological tissues. A theoretical investigation also demonstrates the significant superiority of our approach compared to confocal microscopy or OCT. In particular, an imagingdepth limit of 22 is predicted for smartOCT in biological tissues, hence pushing back drastically the fundamental multiple scattering limit in optical imaging.
Materials and Methods
The following components were used in the experimental setup [Fig. 1]: A femtosecond laser (FEMTOSECOND^{TM} FUSION^{TM} 20400), a spatial light modulator (PLUTO NIR2, Holoeye), an objective lens (Olympus, 10, NA=0.25), a CCD camera (DALSA Pantera 1M60) with a dynamic range of 60 dB. The incident light power is of 10 mW in the experiment. The radiant flux is thus of W.cm at the focal spot in free space. For each wavefront input, the complex reflected wavefield is extracted from four intensity measurements. Hence, the measurement of each line of the reflection matrix can be done at 15 Hz. According to the singletomultiple scattering ratio, the reflected wavefield has to be averaged over a given number of measurements. In the imaging experiments through the paper layer [Figs. 2 and 3] and biological tissues [Fig. 4], this number has been fixed to 32 and 5, respectively. Hence, the duration time for the recording of the reflection matrix has been of 10 min for the paper (289 input wavefronts) and 5 min for biological tissues (961 input wavefronts). The numerical postprocessing of the reflection matrix (single scattering projection and iterative timereversal) to get the final image only takes a few seconds.
Acknowledgments
Fundings The authors are grateful for funding provided by LABEX WIFI (Laboratory of Excellence within the French Program Investments for the Future, ANR10LABX24 and ANR10IDEX000102 PSL*). A. B. acknowledges financial support from the French “Direction Générale de l’Armement”(DGA). D. L. acknowledges financial support from LABEX WIFI and European Research Council (ERC Synergy HELMHOLTZ). G. L. and A. A. would like to acknowledge funding from High Council for Scientific and Technological Cooperation between France and Israel under reference P2R Israel N 29704SC.
Author Contributions A.A. initiated and supervised the project. A.A. conceived the experiment. A.B. and D.L. built the experimental set up and performed the experiments. A.B., D.L. and A.A. analyzed the experiments. A.A. performed the theoretical study. A.B. and A.A. prepared the manuscript. All authors discussed the results and contributed to finalizing the manuscript.
Supplementary Materials
.1 Optical characterization of the paper layers
White paper sheets are used as turbid media in our experiment. White paper is a strongly scattering medium with negligible absorption Hubbe (). In this section, we report on the measurement of the transport parameters that govern the diffusion of light across each paper sheet.
Scattering mean free path
The ballistic component of the total transmitted intensity decays exponentially across a scattering layer of thickness akkermans (),
(S1) 
with the incident intensity. The scattering mean free path can be measured by investigating the attenuation of the ballistic light across an increasing number of paper sheets. To that aim, the experiment described in Fig. S1A has been performed. A collimated laser beam (=810 nm) illuminates the paper sheets. The transmitted intensity is the sum of the ballistic and the diffuse light. A pinhole is placed at the output of the medium to spatially discriminate the ballistic photons. The transmitted intensity is recorded on a CCD camera placed 50 cm behind the pinhole.
The transmitted intensity pattern is shown in Fig. S1B for configurations with 1, 2 and 3 sheets of paper. The left column represents a typical intensity pattern recorded for one realization of disorder. The right column represents the average of the transmitted intensity over a given number of disorder realizations. This is done by scanning spatially the surface of the scattering medium. For 1 and 2 layers, the ballistic component emerges on top of the diffusive halo. Its intensity can be estimated by subtracting the maximum intensity with the background intensity in the vicinity of the ballistic peak. For 3 layers of paper, the ballistic component cannot be revealed despite averaging over realizations of disorder. This confirms that white paper is strongly scattering. The thickness of each paper sheet has been measured with a precision caliper: m. The optical thickness of one paper sheet can be obtained by fitting linearly as a function of the number of layers [Fig. S1C]. We obtain . This leads to the following estimation of the scattering mean free path: 13.4 m.
Transport Mean Free Path
Since , most of the transmitted light is scattered several times while propagating through one sheet. For several sheets, the diffusion approximation becomes valid. The total transmission coefficient can then be written as follows durian (),
(S2) 
with the distance from the scattering layer boundary at which the energy flux should cancel according to diffusion theory. In presence of a refractive index mismatch at the scattering madium boundaries, it can be expressed as Xhu ()
(S3) 
with the internal reflection coefficient. Considering a mean refractive index of 1.5 for the paper Hubbe (), Xhu () and .
The measurement of has been performed with the experiment described in Fig. S2A. A collimated laser beam illuminates the scattering sample and most of the transmitted light is collected with a powermeter. The transmission coefficient is measured for an increasing number of paper layers. Figure S2B displays as a function of M. A linear fit of experimental points leads to an estimation of using Eq.S2. With , we find 19.9m.
.2 Singletomultiple scattering ratio for various imaging techniques
To address the issue of MS in optical imaging, it is important to establish theoretically the limit of existing imaging techniques in inhomogeneous media. To that aim, the relevant parameter is the singletomultiple scattering ratio (SMR). In the following, this quantity is derived theoretically for several imaging techniques (conventional/confocal microscopy, OCT, smartOCT) to highlight the gain in SMR provided by each of them. To that aim, we will consider the detection and imaging of an object embedded in or hidden behind a strongly scattering medium at a depth .
Let us first consider a conventional microscopy configuration. The target is placed in the focal plane of a microscope objective (MO). Under a plane wave illumination, the SMR of the backscattered wavefield in the backfocal plane of a MO is given by,
(S4) 
with , the differential scattering crosssection of the target, , the fieldofview (FOV) of the optical system and , the static albedo of the scattering layer akkermans (); Akkermans2 (). Not surprisingly, SS is favored by the brightness of the target. However, the most important parameter here is the attenuation undergone by the ballistic wave across the scattering layer. Only a very tiny fraction of the incident energy, , is converted into the SS contribution used for imaging. In addition to the severe attenuation of the ballistic wave, the turbid medium gives rise to a speckle wavefield whose intensity is given by the static albedo . For a scattering layer of thickness , its expression is given by akkermans ()
(S5) 
When one tries to image the target with conventional microscopy, the SMR in the conjugate image plane, referred to as SMR, can be expressed as follows
(S6)  
(S7) 
with , the Strehl ratio ranging from 0 to 1 Mahajan (), , the resolution length of the imaging system and , the MO numerical aperture. Compared to conventional microscopy, the SMR is increased by a factor which corresponds the number of resolution cell in the FOV. Whereas the MS background results from the incoherent superimposition of independent speckle grains, the target image results from the constructive interference of the ballistic photons over the numerical aperture. However, the aberration undergone by the target ballistic wavefront across the scattering layer degrades the target focal spot and lowers its intensity. This is accounted for by the Strehl ratio in Eq.S6. is directly proportional to the focusing parameter introduced by Mallard and Fink in the ultrasound imaging context mallart ().
To cope with the fundamental issue of MS, several approaches have been proposed in order to enhance the SS contribution drowned into a predominant MS background. The first option is to spatially discriminate SS and MS as performed in confocal microscopy. Ideally the incoming radiation is focused to a single voxel and only light backscattering from that voxel is collected, allowing to reject a large number of multiplyscattered photons . However, scattered light can blur the focused beam outside the target volume and unwanted photons from other voxels can be scattered back into trajectories that will be collected by the microscope. Theoretically, the SMR provided by confocal microscopy, referred to as SMR, can be expressed as,
(S8)  
(S9) 
Compared to conventional microscopy, the SMR is increased a new time by the factor which results from the coherent summation of the incident ballistic wavefront at the target location. The Strehl ratio in Eq.S8 accounts for the aberration effect undergone by this incident wavefront.
The second way to enhance SS relatively to MS consists in discriminating SS from MS photons with coherence time gating. It allows to select the ballistic photons over a time window centered around their time of flight. Probably, the most widely employed coherence timegated technique is optical coherence tomography (OCT), which combines scanning confocal microscopy with coherent heterodyne detection. The MS intensity is now given by a timedependent albedo akkermans (); Akkermans2 (), with , the ballistic time. The SMR in OCT can be deduced from Eq.S9 by substituting the static albedo with the ratio , with the bandwidth of the light source (inversely proportional to its coherence time ). The SMR in OCT, referred to as SMR, is thus given by,
(S10)  
(S11) 
Compared to confocal microscopy, the SMR in OCT is increased by a factor which accounts for the number of multiply scattered photons rejected by coherence time gating. The timedependent albedo differs according to the imaging configuration. If the target is placed behind a scattering layer of thickness , the following expression of should be considered Lagendijk ()
(S12) 
where is the diffusion constant that governs wave transport across the scattering layer. If the target is embedded within the scattering medium, the timedependent albedo for a semiinfinite medium should be considered akkermans (); Akkermans2 ()
(S13) 
According to the experimental configuration, the timedependent albedo either displays a power law or an exponential decrease with the timeofflight. In both cases, coherence time gating allows to drastically reject multiplyscaterred photons.
In smartOCT, in addition to confocal and coherence time gating operations, an eigenvalue decomposition of the single scattering reflection matrix allows to get rid of the residual multiple scattering contribution. Whereas each ballistic echo emerges along one single eigenstate of (the signal subspace), the incoherent MS wavefield emerges with the same probability along its eigenstates aubry2 (). This implies an enhancement of the SMR by a factor compared to OCT [Eq.S10]. Moreover, smartOCT does not suffer from aberration issues in terms of detection. Timereversal processing directly yields the distorted wavefront that compensates for the aberration effects induced by the scattering layer prada2 (); popoff3 (); aubry2 (). The SMR ratio is thus independent of the Strehl ratio . As a consequence, the SMR associated with smartOCT, referred to as SMR, can be expressed as,
(S14)  
This theoretical study is first applied to the experimental configuration described in the accompanying paper. The parameters used for the computation of the SMR are described in Tab.S1. The transport parameters of the paper (, , ) have been derived in the first section of this supplementary material. The experimental parameters (, , NA, , ) correspond to those reported in the paper. The differential scattering crosssection has been estimated from Mie theory bohren1983absorption () by considering a 10 mdiameter ZnO spherical bead as a target. The numerical value given in Tab.S1 corresponds to an average of over the numerical aperture of the MO. At last, the Strehl ratio is estimated from , the coherence length of the ballistic wavefield in the focal plane, such that . As m in the reported experiment and m, it yields .
Figure S3 displays the evolution of the SMR as a function of the optical thickness by applying numerically Eqs.S4, S6, S8, S10 and S14 with the parameters shown in Tab.S1. Note that Eq.S12 is considered for the computation of the timedependent albedo . Figure S3 illustrates how a confocal illumination and a coherence time gating allows to drastically improve the SMR compared to conventional microscopy. Nevertheless, the theoretical OCT imagingdepth limit (SMR) remains limited to in this configuration. This explains why the timegated confocal image is not able to reveal the presence of a target for an optical thickness 12 in our experiment. On the contrary, the smartOCT imagingdepth limit is predicted to be around 12 . This is in a remarkable agreement with our experimental results showing a successful target detection for an optical tickness . Contrastedly, the target is far to be detectable and imaged whether it be by conventional microscopy (SMR , Eq.S6), confocal microscopy (SMR , Eq.S8) or by OCT (SMR , Eq.S10). The SMR of the backscattered wavefield is of (Eq.S4) meaning that only 1 reflected photon over one thousand billions is associated to a single scattering event from the target in our experiment.
In Fig. 5, the SMR is displayed as a function of the optical depth in the context of biological tissues imaging. Note that the SMR for conventional and confocal microscopy (Eqs.S6S8) has been computed by considering the asymptotic limit of the static albedo (, see Eq.S5). As to coherence time gating, Eq.S13 is considered for the computation of the timedependent albedo as the target is assumed to be embedded within the scattering medium. The parameters used for the computation of the SMR are described in Tab.S1. The considered transport parameters (, ) are typical of invivo cortex tissues Schott (). The experimental parameters (, , NA, , ) are typical of fullfield OCT Assayag (). The internal reflection coefficient is assumed to be zero since the use of an immersion microscope objective will limit the impedance mismatch with tissues. The target scattering crosssection is arbitrarily chosen to be the same as in the reported experiment. At last, the Strehl ratio in brain tissues is estimated from a two photons microscopy experiment that reports a fivefold signal enhancement when optical aberrations from the brain tissues are corrected with adaptive optics Ji ().
Experimental  Reported experiment:  Fullfield oct: 

configuration  paper sheets Hubbe ()  Brain tissues Assayag () 
1.5 Hubbe ()  1.4 Jacques ()  
[m ]  13.4  200 Schott () 
[m ]  19.9  2000 Schott () 
0.57 Xhu ()  0  
[m.sr]  
0.1  0.4  
[mm]  1  na 
[nm]  810  810 
NA  0.25  0.4 
[m]  40  1000 
[fs]  50  5 
.3 Eigenvalue distribution of the reflection matrix in the multiple scattering regime
In this section, we derive the numerical method to obtain the eigenvalue distribution of expected in a fully multiple scattering regime [see Fig. 2F]. If the coefficients of were complex random variables independently and identically distributed, the eigenvalue distribution would follow the MarcenkoPastur law marcenko (); tulino (). However, this assumption is not fulfilled here. First, all the elements of do not exhibit the same variance because of the single scattering filtering operation [see Figs. 2E and 3E]. Second, some residual correlations may arise in the measured matrix due to experimental imperfections. In particular, a slight curvature of the reference beam in the experiment may induce some shortrange correlations of the wavefield at the output.
The correlation between two coefficients and of can be expressed as sengupta ()
(S15) 
where the symbol denotes an ensemble average. and are matrices. We will refer to them as the correlation matrices. As the correlation properties are statistically invariant by translation, and are Toeplitz matrices: and . The correlation coefficients, and , between the columns and the lines of can be estimated as follows
(S16) 
where the symbol denotes an average over the variables in the subscript. From these correlation coefficients, one can estimate the correlation matrices and .
Once and are estimated, one can deduce the eigenvalue distribution expected in a fully multiple scattering regime aubry_wrmc (). It consists in generating numerically a matrix whose elements are circularly symmetric complex gaussian random variables with zero mean. Then, a matrix is built from , such that
(S17) 
One can show that the matrix exhibits the same correlation properties at emission and reception as the experimental matrix . The matrix is then projected on the characteristic SS matrix [Eq.3]
(S18) 
A singular value decomposition of is then performed. The eigenvalues of (square of the singular values) are renormalized by their mean. An histogram of the eigenvalues is then obtained by averaging over 2000 numerically generated random matrices . The resulting theoretical distribution is shown in Fig. 2F.
References
 A. P. Mosk, A. Lagendijk, G. Lerosey, M. Fink, Controlling waves in space and time for imaging and focusing in complex media. Nature Photon. 6, 283292 (2012).
 A. Derode, P. Roux, M. Fink, Robust acoustic time reversal with highorder multiple scattering. Phys. Rev. Lett. 75, 42064209 (1995).
 A. Derode, A. Tourin, J. de Rosny, M. Tanter, S. Yon, M. Fink, Taking advantage of multiple scattering to communicate with timereversal antennas. Phys. Rev. Lett. 90, 014301 (2003).
 I. M. Vellekoop, A. P. Mosk, Focusing coherent light through opaque strongly scattering media. Opt. Lett. 32, 2309–2311 (2007).
 S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, 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).
 M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.H. Park, W. Choi, Maximal energy transport through disordered media with the implementation of transmission eigenchannels. Nature Photon. 6, 583587 (2012).
 I. N. Papadopoulos, S. Farahi, C. Moser, D. Psaltis, Focusing and scanning light through a multimode optical fiber using digital phase conjugation. Opt. Exp. 20, 1058310590 (2012).
 S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, S. Gigan, Image transmission through an opaque material. Nat. Commun. 1, 1–5 (2010).
 T. Cizmar, K. Dholakia, Exploiting multimode waveguides for pure fibrebased imaging. Nature Comm. 3, 1027 (2012).
 T. Durduran, R. Choe, W. B. Baker, A. G. Yodh, Diffuse optics for tissue monitoring and tomography. Rep. Prog. Phys. 73, 076701 (2010).
 S. G. Resink, A. C. Boccara, W. Steenbergen, Stateofthe art of acoustooptic sensing and imaging of turbid media. J. Biomed. Opt. 17, 040901 (2012).
 M. Xu, L. V. Wang, Photoacoustic imaging in biomedicine. Rev. Sci. Instrum. 77, 041101 (2006).
 I. Freund, M. Rosenbluh, S. Feng, Memory effects in propagation of optical waves through disordered media. Phys. Rev. Lett. 61, 23282331 (1988).
 S. Feng, C. Kane, P. A. Lee, A. D. Stone, Correlations and fluctuations of coherent wave transmission through disordered media. Phys. Rev. Lett. 61, 834837 (1988).
 J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, A. P. Mosk, Noninvasive imaging through opaque scattering layers. Nature 491, 232234 (2012).
 O. Katz, E. Smal, Y. Silberberg, Looking around corners and through thin turbid layers in real time with scattered incoherent light. Nature Photon. 6, 549553 (2012).
 O. Katz, P. Heidmann, M. Fink, S. Gigan, Noninvasive singleshot imaging through scattering layers and around corners via speckle correlations. Nature Photon. 8, 784790 (2014).
 J. B. Pawley, ed., Handbook Of Biological Confocal Microscopy (SpringerVerlag, New York, 2006).
 V. Ntziachristos, Going deeper than microscopy: the optical imaging frontier in biology. Nature Meth. 7, 603614 (2010).
 P. Theer, M. T. Hasan, W. Denk, Twophoton imaging to a depth of 1000 m in living brains by use of a ti:al2o3 regenerative amplifier. Opt. Lett. 28, 10221024 (2003).
 J. G. Fujimoto, S. De Silvestri, E. P. Ippen, C. A. Puliafito, R. Margolis, A. Oseroff, Femtosecond optical ranging in biological systems. Opt. Lett. 11, 150152 (1986).
 L. Wang, P. P. Ho, C. Liu, G. Zhang, R. R. Alfano, Ballistic 2d imaging through scattering walls using an ultrafast optical kerr gate. Science 253, 769771 (1991).
 G. E. Anderson, F. Liu, R. R. Alfano, Microscope imaging through highly scattering media. Opt. Lett. 19, 981983 (1994).
 D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, J. G. Fujimoto, Optical coherence tomography. Science 254, 11781181 (1991).
 J. G. Fujimoto, M. E. Brezinski, G. J. Tearney, S. A. Boppart, B. Bouma, M. R. Hee, J. F. Southern, E. A. Swanson, Optical biopsy and imaging using optical coherence tomography. Nature Med. 1, 970972 (1995).
 A. F. Fercher, Optical coherence tomography. J. Biomed. Opt. 1, 157173 (1996).
 W. Drexler, J. G. Fujimoto, eds., Optical Coherence Tomography. Technology and Applications. (SpringerVerlag, New York, 2008).
 P. Andersen, T. M. Joørgensen, L. Thrane, A. Tycho, H. T. Yura, Optical coherence tomography. Technology and Applications (SpringerVerlag, New York, 2008, 2008), chap. Modeling LightTissue Interaction in Optical Coherence Tomography Systems, p. 73.
 C. Dunsby, P. W. French, Techniques for depthresolved imaging through turbid media including coherencegated imaging. J. Phys. D: Appl. Phys. 36, R207R227 (2003).
 A. Aubry, A. Derode, Random Matrix Theory Applied to Acoustic Backscattering and Imaging In Complex Media. Phys. Rev. Lett. 102, 084301 (2009).
 A. Aubry, A. Derode, Detection and imaging in a random medium : A matrix method to overcome multiple scattering and aberration. J. Appl. Phys. 106, 044903 (2009).
 S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, A. Derode, A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering. Appl. Phys. Lett. 104, 234105 (2014).
 Y. Choi, T. R. Hillman, W. Choi, N. Lue, R. R. Dasari, P. T. C. So, W. Choi, Z. Yaqoob, Measurement of the timeresolved reflection matrix for enhancing light energy delivery into a scattering medium. Phys. Rev. Lett. 111, 243901 (2013).
 S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J.S. Lee, Y.S. Lim, Q. H. Park, W. Choi, Imaging deep within a scattering medium using collective accumulation of singlescattered waves. Nature Photon. 9, 253258 (2015).
 J.L. Robert, 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, 866877 (2008).
 J.L. Robert, M. Fink, The timereversal operator with virtual transducers: Application to farfield aberration correction. J. Acoust. Soc. Am. 124, 36593668 (2008).
 C. Prada, M. Fink, Eigenmodes of the time reversal operator: A solution to selective focusing in multipletarget media. Wave motion 20, 151–163 (1994).
 C. Prada, S. Manneville, D. Spoliansky, M. Fink, Decomposition of the time reversal operator: Detection and selective focusing on two scatterers. J. Acoust. Soc. Am. 99, 20672076 (1996).
 S. M. Popoff, A. Aubry, G. Lerosey, M. Fink, A. C. Boccara, S. Gigan, Exploiting the timereversal operator for adaptive optics, selective focusing, and scattering pattern analysis. Phys. Rev. Lett. 107, 263901 (2011).
 D. Psaltis, A. S. Goy, Digital confocal microscope. Opt. Exp. 20, 2272022727 (2012).
 J.L. Robert, M. Fink, The prolate spheroidal wave functions as invariants of the time reversal operator for an extended scatterer in the fraunhofer approximation. J. Acoust. Soc. Am. 125, 218226 (2008).
 O. Assayag, M. Antoine, B. SigalZafrani, M. Riben, F. Harms, A. Burcheri, B. Le Conte de Poly, A. C. Boccara, Large field, high resolution full field optical coherence tomography: A preclinical study of human breast tissue and cancer assessment. Technol Cancer Res Treat. 13, 455468 (2014).
 A. Badon, G. Lerosey, A. C. Boccara, M. Fink, A. Aubry, Retrieving timedependent green’s functions in optics with lowcoherence interferometry. Phys. Rev. Lett. 114, 023901 (2015).
 M. A. Hubbe, J. J. Pawlak, K. A.A., Paper’s appearance: A review. BioRes. 3, 627665 (2008).
 E. Akkermans, G. Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University Press, London, 2006).
 D. J. Durian, Influence of boundary reflection and refraction on diffusive photon transport. Phys. Rev. E 50, 857866 (1994).
 J. X. Xhu, D. J. Pine, D. A. Weitz, Internal reflection of diffusive light in random media. Phys. Rev. A 44, 39483959 (1991).
 E. Akkermans, P. E. Wolf, R. Maynard, G. Maret, Theoretical study of the coherent backscattering of light by disordered media. J. Phys. France 49, 7798 (1988).
 V. N. Mahajan, Strehl ratio for primary aberrations: some analytical results for circular and annular pupils. J. Opt. Soc. Am. 72, 12581266 (1982).
 R. Mallart, M. Fink, Adaptive focusing in scattering media through soundspeed inhomogeneities: The van cittert zernike approach and focusing criterion. J. ACoust. Soc. Am. 96, 37213732 (1994).
 A. Lagendijk, R. Vreeker, P. de Vries, Influence of internal reflection on diffusive transport in strongly scattering media. Phys. Lett. A 136, 8188 (1989).
 C. F. Bohren, D. R. Huffman, Absorption and scattering of light by small particles (John Wiley & Sons, Inc., New York, 1983).
 S. Schott, J. Bertolotti, J.F. Léger, L. Bourdieu, S. Gigan, Characterization of the angular memory effect of scattered light in biological tissues. Opt. Exp. 23, 1350513516 (2015).
 N. Ji, T. R. Sato, E. Betzig, Characterization and adaptive optical correction of aberrations during in vivo imaging in the mouse cortex. Proc. Natl. Acad. Sci. USA 109, 2227 (2012).
 M. Marc̆enko, L. Pastur, Distributions of eigenvalues for some sets of random matrices. Math. USSRSbornik 1, 457483 (1967).
 A. Tulino, S. Verdù, Random matrix theory and wireless communications. Fundations and Trends in Communications and Information Theory 1, 1182 (2004).
 A. M. Sengupta, P. P. Mitra, Distribution of singular values for some random matrices. Phys. Rev. E 60, 33893392 (1999).
 A. Aubry, A. Derode, Singular value distribution of the propagation matrix in random scattering media. Waves Random Complex Media 20, 333363 (2010).
 S. L. Jacques, Optical properties of biological tissues: a review. Phys. Med. Biol. 58, R37R61 (2013).