Deep imaging inside scattering media through virtual spatiotemporal wavefront shaping

Light scattering from inhomogeneities scrambles the wavefront and attenuates the signal exponentially while obscuring it with a speckled multiple-scattering background ^{1, 2}. Overcoming light scattering to image at high resolution deep inside opaque media is a grand challenge that can enable through-skull optical imaging and stimulation, surgery-free biopsy, observation of cellular processes in their natural state, optical characterization of bulk colloids, high-resolution nondestructive device testing, and many other applications.

Traditional imaging methods tackle this problem by filtering away the multiple-scattering background. The reflectance confocal microscopy (RCM) ^{3, 4} uses spatial gating. Optical coherence tomography (OCT) ⁵ applies a coherence gate. The axial span of the spatial beam (which is twice the Rayleigh range) must cover the axial span of the image (i.e., the depth of field $\Delta z$), so OCT operates with a small numerical aperture (NA), leading to a coarse lateral resolution (typically 5–20 µm). Optical coherence microscopy (OCM) ⁶ combines spatial filtering and coherence gating while trading off the depth of field. Interferometric synthetic aperture microscopy ⁷ extends the depth of field, though its signal decreases with the loss of confocal gates away from the focal plane. Hardware ⁸, computational ^{9, 10, 11, 12}, and hybrid ¹³ adaptive optics can restore the resolution of images blurred by aberrations but cannot overcome severe scattering as they do not control the numerous scattering modes in both the incident and the return paths. Multiphoton microscopy adopts longer wavelengths where scattering is weaker but requires fluorescence labeling and high-power lasers ^{14, 15}. Tomographic phase microscopy quantitatively obtains the refractive index distribution ^{16, 17, 18, 19, 20} but requires thin specimens and access to both sides of the sample for transmission measurements. Selective detection of the forward scattered photons ^{21, 22} and photoacoustic tomography ²³ allow imaging deeper at reduced resolution. The depth-over-resolution ratio stays around or below 200 in biological tissue across most label-free imaging methods (including non-optical ones such as ultrasound, X-ray CT, and MRI) ^{23, 24, 25, 26}.

Over the past decade, a paradigm shift came with the appreciation that light scattering is deterministic and controllable. Given enough control elements, a suitable incident wavefront can reverse multiple scattering ^{27, 28} thanks to the time-reversal symmetry of Maxwell’s equations. Adopting such wavefronts can, in principle, allow imaging at much deeper depths. But the initial methods to determine these special wavefronts come with serious limitations as they introduce physical guidestars at the imaging site ²⁹ such as invasive cameras ^{30, 31} or photorefractive crystals ^{32, 33}, fluorescence labels ^{34, 35, 36, 37}, acoustic foci with greatly reduced spatial resolution ^{38, 39}, or isolated moving particles ^{40, 41}. Additionally, each wavefront is only valid within a volume called the isoplanatic patch ⁸, whose size shrinks with depth and is less than 10 µm at 1 mm deep inside biological tissue ^{42, 43}. Most of these methods also require a spatial light modulator (SLM) with a limited update speed. Such physical-guidestar-based schemes are schematically illustrated in Fig. 1a. To avoid using physical guidestars, researchers developed incoherent imaging methods based on angular correlations ^{44, 45} and optimization ^{46, 47, 48, 49}. However, they only apply to planar targets at a large distance away from the scattering media because there is no gating mechanism for depth sectioning and because the small angular memory effect range limits the field of view. Reflection matrix measurement is another way to avoid physical guidestars. One can accumulate single-scattering signals using correlations in the reflection matrix ⁵⁰ and remove aberrations and scattering through singular value decomposition ^{51, 52} or iterative phase conjugation ^{53, 54, 55, 56, 57, 58}, culminating in the “volumetric reflection-matrix microscopy” (VRM) method ⁵⁸ that performs both spectral and angular corrections at multiple depths. These matrix-based methods show great potential but, to date, have yet to deliver the drastic beyond-adaptive-optics scattering reversal promised by wavefront shaping.

Here we use the scattering matrix to unify different approaches and combine the advantages of both the traditional and the recent methods, realizing label-free volumetric high-resolution imaging across a large depth of field deep inside the scattering media, with multi-isoplanatic spatiotemporal wavefront corrections on top of confocal and temporal gates while avoiding the restrictions of physical guidestars and slow SLMs. As illustrated in Fig. 1b, we measure the hyperspectral scattering matrix of the sample and use it to synthesize a spatiotemporal focus with pulse compression, refractive index mismatch correction, and wavefront correction for every isoplanatic volume. The focus is scanned digitally to yield a 3D image of the sample with no trade-off between the lateral resolution and the depth of field. We use the reconstructed image as virtual guidestars intrinsic to the sample and develop an optimization with regularization and progression strategies to find the optimal dispersion compensation and input/output wavefronts. This enables a depth-over-resolution ratio of 910 when imaging a resolution target at one millimeter beneath ex vivo mouse brain tissue—the highest reported in the literature (including both optical and non-optical label-free methods) to our knowledge. We attain diffraction-limited resolution when the signal is reduced by over ten-million-fold due to multiple scattering. We maintain the ideal transverse and axial resolutions across a depth of field of over 70 times the Rayleigh range at depths beyond three transport mean free paths inside an opaque colloid—the deepest reported to our knowledge. By formulating the problems of imaging and reversing light scattering as reconstruction and optimization problems, this approach offers a versatile platform that uses wave physics and computation to reveal hidden information.

The scattering matrix $S(\bf k_{\rm out},\bf k_{\rm in},\omega)$ encapsulates the sample’s complete linear response ^{59, 60} (Fig. 1c). Any incident wave is a superposition of plane waves: ${\it E}_{\rm in}({\bf r},t)=\sum_{{\bf k}_{\rm in},\omega}\tilde{\it E}_{\rm in}({\bf k}_{\rm in},\omega)e^{i{\bf k}_{\rm in}\cdot{\bf r}-i\omega t}$, where ${\bf r}=(x,y,z)$ is the position, $t$ is time, and $\tilde{\it E}_{\rm in}({\bf k}_{\rm in},\omega)$ is the amplitude of its plane-wave component with incoming momentum $\bf k_{\rm in}$ at frequency $\omega$. The amplitudes $\tilde{E}_{\rm out}({\bf k}_{\rm out},\omega)$ that form the resulting outgoing wave ${\it E}_{\rm out}({\bf r},t)=\sum_{{\bf k}_{\rm out},\omega}\tilde{\it E}_{\rm out}({\bf k}_{\rm out},\omega)e^{i{\bf k}_{\rm out}\cdot{\bf r}-i\omega t}$ are given by the scattering matrix through $\tilde{E}_{\rm out}({\bf k}_{\rm out},\omega)=\sum_{{\bf k}_{\rm in}}S({\bf k}_{\rm out},{\bf k}_{\rm in},\omega)\tilde{E}_{\rm in}(\bf k_{\rm in},\omega)$. The angular summations are restricted to $|{\bf k}|^{2}=|{\bf k}_{\parallel}|^{2}+k_{z}^{2}=(n_{\rm bg}\omega/c)^{2}$ in a background medium with speed of light $c/n_{\rm bg}$; we use “angle” interchangeably with “momentum.”

After measuring a subset of $S(\bf k_{\rm out},\bf k_{\rm in},\omega)$, one can digitally synthesize the sample’s response with tailored measurements in space-time given customized spatiotemporal inputs, which we collectively refer to as “virtual spatiotemporal wavefront shaping.” As illustrated in Fig. 1d, summing over incident angles $\bf k_{\rm in}$ with incoming plane-wave amplitudes $\tilde{E}_{\rm in}({\bf k}_{\rm in},\omega)=e^{-i{\bf k_{\rm in}}\cdot{\bf r_{\rm in}}}$ creates an incident beam ${\it E}_{\rm in}({\bf r},\omega)$ spatially focused at $\bf r_{\rm in}$, forming an input spatial gate; summing over outgoing angles $\bf k_{\rm out}$ with outgoing plane-wave amplitudes $\tilde{E}_{\rm out}({\bf k}_{\rm out},\omega)$ from the scattering matrix yields the scattered field ${\it E}_{\rm out}({\bf r}={\bf r}_{\rm out},\omega)$ given a virtual detector at position ${\bf r}_{\rm out}$, forming an output spatial gate; summing over frequency $\omega$ with an $e^{-i\omega(t-t_{0})}$ phase creates an incident pulse that arrives at $\bf r_{\rm in}$ at time $t=t_{0}$, forming a temporal gate. Given a spatiotemporal focus arriving at $\bf r_{\rm in}$ at time $t_{0}$, the scattered field is therefore $E_{\rm out}({\bf r}_{\rm out},\it t)=\sum_{\bf k_{\rm out},\bf k_{\rm in},\omega}S({\bf k}_{\rm out},{\bf k}_{\rm in},\omega)e^{i{\bf k}_{\rm out}\cdot{\bf r}_{\rm out}-i{\bf k}_{\rm in}\cdot{\bf r}_{\rm in}-i\omega(t-t_{0})}$. We then perform a virtual experiment with a confocal spatial focus (${\bf r}_{\rm in}={\bf r}_{\rm out}={\bf r}$) and with the temporal focus aligned with the spatial one (evaluating $E_{\rm out}$ at time $t=t_{0}$), which yields the triply-gated scattering amplitude of the sample at position ${\bf r}$, denoted as

Digitally scanning $\bf{r}$ (Fig. 1b) forms a phase-resolved image of the sample where the input spatial gate, output spatial gate, and temporal gate all align at every point in the absence of aberrations and scattering. The lateral spread of $\bf k_{\rm out}-\bf k_{\rm in}$ in the summation (typically set by the NA) determines the lateral resolution, with no degradation away from any particular focal plane. The axial spread (typically set by the spectral bandwidth) determines the axial resolution. The triple summation averages away random noises, providing a signal-to-noise ratio advantage akin to that of frequency-domain OCT over time-domain OCT ⁶¹. We use the non-uniform fast Fourier transform (NUFFT) ⁶² to efficiently evaluate these summations and the 3D spatial scan. Eq. (1) is the minimal, optimization-free version of what we call “scattering matrix tomography” (SMT).

VRM used an object momentum conservation principle ^{50, 58} to reconstruct images at multiple en-face slices by performing time gating at one depth, summing over incident angles while fixing the in-plane momentum transfer, performing a 2D inverse Fourier transform, and then looping over the depth ⁵⁸. In comparison, SMT uses a virtual scanning of a confocal spatiotemporal focus, which results in a simpler and more intuitieve formalism, enables slicing in any orientation or direct volumetric reconstruction, generalizes the use of reflection matrix to any subset of the scattering matrix (can also be transmission, remission ^{21, 63}, or a combination of them), adopts NUFFT which is much faster than brute-force summations and loops, and offers a rigorous basis for performing the index-mismatch correction and focus-quality-based spatiotemporal wavefront corrections below.

The triple gating suppresses multiple scattering. The single-scattering field at output ${\bf k}_{\rm out}$ for an input from ${\bf k}_{\rm in}$ is given by the Born approximation as proportional to $\tilde{\eta}({\bf q})$, namely the ${\bf q}={\bf k}_{\rm out}-{\bf k}_{\rm in}$ component of the 3D Fourier transform of the sample’s permittivity contrast profile $\eta({\bf r})$ ^{64, 16, 17, 5}. Such single-scattering contributions add up in phase in Eq. (1) to form the image, similar to an inverse Fourier transform from $\tilde{\eta}({\bf q})$ to $\eta({\bf r})$. The multiple-scattering contributions do not add up in phase. Therefore, the triple summations over $\bf k_{\rm out}$, $\bf k_{\rm in}$, and $\omega$ boost the single-to-multiple-scattering ratio in $\psi_{\rm SMT}({\bf r})$ compared to the raw measurement data $S({\bf k}_{\rm out},{\bf k}_{\rm in},\omega)$.

To remove aberrations and scattering, we perform digital spatiotemporal wavefront corrections. The chromatic aberrations in the optical elements of the system and the frequency dependence of the refractive index in the sample create a frequency-dependent phase that misaligns and broadens the temporal gate. In SMT, we digitally compensate for such dispersion using a spectral phase $\theta(\omega)$ acting as a virtual pulse shaper (Fig. 1e), similar to the dispersion compensation in spectral-domain OCT ⁶⁵. The refractive index mismatch between the sample medium (e.g., biological tissue), the far field (e.g., air), and the coverslip (if there is one) refracts the rays, degrades the input and output spatial gates ⁶⁶, and misaligns the spatial gates and the temporal gate. By using the appropriate propagation phase shift $(\bf k_{\rm out}-\bf k_{\rm in})\cdot{\bf r}$ in each medium and the appropriate arrival time, SMT achieves an ideal spatiotemporal focus at all depths even in the presence of refraction, effectively creating a virtual dry objective lens that reverses both the spatial and the chromatic aberrations arising from the index mismatch (Fig. 1f) without liquid immersion (Supplementary Sect. III C). Importantly, we also digitally introduce angle-dependent phase profiles $\phi_{\rm in}(\bf k_{\rm in})$ and $\phi_{\rm out}(\bf k_{\rm out})$ that act as two virtual SLMs, one for the incident wave and one for the outgoing wave (Fig. 1g), which can adopt any wavefront to compensate for scattering and additional aberrations. This yields the general form of SMT,

By measuring the reflection from a mirror, we remove the dispersion and aberrations in the input path of the optical system (Supplementary Sect. III B).

To find the digital corrections $\theta(\omega)$, $\phi_{\rm in}(\bf k_{\rm in})$, and $\phi_{\rm out}(\bf k_{\rm out})$, we recognize that a good spatiotemporal focus at $\mathbf{r}$ enhances the scattering intensity from targets at that position, so $I_{\rm SMT}(\bf{r})$ acts as a virtual noninvasive guidestar at $\mathbf{r}$ (Fig. 1b). Therefore, we maximize an image quality metric $M$,

with $I_{0}$ being a constant; the logarithm factor promotes image sharpness ^{67, 12}. Image-metric-based optimizations have long been used in adaptive optics ^{8, 9, 10, 11, 12, 13} and more recently with wavefront shaping without a time gate ^{46, 48, 37}; motivated by the versatility of optimization problems, here we generalize this approach to perform double-path spatiotemporal digital wavefront correction. We derive the gradient of $M$ with respect to the parameters and adopt a quasi-Newton method, the low-storage Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method ⁶⁸, in the NLopt library ⁶⁹. Since a wavefront $\phi_{\rm in}(\bf k_{\rm in})$ or $\phi_{\rm out}(\bf k_{\rm out})$ can only remain optimal within an isoplanatic volume ^{8, 42, 43}, we divide the space into zones (both axially and laterally) and use different $\phi_{\rm in/out}$ for different zones.

A direct optimization, however, performs poorly when the scattering is sufficiently strong that signals from the target are initially buried under the speckled multiple-scattering background (Fig. 2a–f). This happens because (1) the problem is high-dimensional and nonconvex with numerous poor local optima, (2) the large number of parameters in $\theta(\omega)$, $\phi_{\rm in}(\bf k_{\rm in})$, and $\phi_{\rm out}(\bf k_{\rm out})$ in the many zones can lead to overfitting, and (3) the optimization may inadvertently enhance unintended contributions like the speckled background instead of signals from a target at the specified position ${\bf r}$. To overcome these challenges, we employ five strategies to regularize the problem, lower its dimension, and guide the search away from poor local optima associated with unintended contributions (Supplementary Sect. III E–F): (i) Restrict $\theta(\omega)$ to a third-order polynomial ^{65, 12}. (ii) Expand $\phi_{\rm in}(\bf k_{\rm in})$ and $\phi_{\rm out}(\bf k_{\rm out})$ in Zernike polynomials ^{70, 10, 11, 12, 13, 8}. (iii) First optimize $\phi_{\rm in/out}$ over the full volume and then progressively optimize over smaller zones while using the previous $\phi_{\rm in/out}$ as the initial guess. (iv) Truncate the scattering matrix to within the respective spatial zone to avoid overfitting. (v) Increase the number of Zernike terms as we progress to smaller zones, since the low-order terms correct for slowly varying aberrations with a large isoplanatic volume, while the high-order terms correct for the sharp variations from localized scattering with a small isoplanatic volume. These strategies work together to enable an effective optimization (Fig. 2g–k).

Since all targets within an isoplanatic volume share the same optimal wavefront while the speckled backgrounds do not, the summation over space ${\bf r}$ in Eq. (3) promotes the target-to-speckle weight in the gradient and makes the optimization landscape smoother. As discussed earlier, the triple gating through summations over $\bf k_{\rm out}$, $\bf k_{\rm in}$, and $\omega$ in Eq. (2) also promotes the target-to-speckle ratio. Without one of these quadruple summations in $\{{\bf r},\bf k_{\rm out},\bf k_{\rm in},\omega\}$ and when scattering is sufficiently strong, the optimization landscape may not reflect the target contribution and may exhibit many local optima associated with the multiple-scattering background (Supplementary Fig. 14).

Previous matrix-based imaging methods used singular value decomposition ⁵² or iterative phase conjugation ^{53, 54, 55, 56, 57, 58} to perform angular wavefront corrections. They were not formulated as an image-metric optimization like Eq. (3) and did not have the above-mentioned regularization and progression strategies, which limited their performance when scattering is strong. The VRM method ⁵⁸ also employed dispersion compensation, but it did so using an iterative phase conjugation scheme with a single angular summation instead of quadruple summations, which further restricted it to situations with less scattering.

We use off-axis holography ⁷¹ to measure the field $\tilde{E}_{\rm out}({\bf k}_{\rm out},\omega)$ scattered from the sample to the incident side (Fig. 3) through a dry objective lens (Mitutoyo M Plan Apo NIR 100X, NA $=$ 0.5). Here, $S({\bf k}_{\rm out},{\bf k}_{\rm in},\omega)$ reduces to the reflection matrix $R({\bf k}_{\rm out},{\bf k}_{\rm in},\omega)$. A CMOS camera (Photron Fastcam Nova S6) captures 64,000 columns of the reflection matrix per second, a dual-axis galvo scanner (ScannerMAX Saturn 5B) scans the incident angle ${\bf k}_{\rm in}$, and a tunable laser (M Squared SolsTiS 1600) scans the frequency $\omega$. The data acquisition here is two to three orders of magnitude faster than previous measurements of broadband scattering matrices ^{60, 72, 58}. The power onto the sample ranges from 0.02 mW to 0.2 mW. We measure around 250 wavelengths $\lambda$ from 740 nm to 940 nm, 2,900 outgoing angles, and 3,900 incident angles within the NA, over a $51\times 51$ µm² area. The detection sensitivity, currently limited by the residual reflection from the objective lens, is 90 dB. See Supplementary Sects. I–II for details.

SMT provides a conceptually and mathematically simple framework that combines the advantages of existing methods while offering versatile capabilities. The main aspect not considered in this work is the frequency dependence of the wavefront correction. It is well known that the optimal wavefront for focusing behind scattering media depends sensitively on the frequency ⁷⁴. Future work can explore optimization strategies that incorporate the frequency dependence in $\phi_{\rm in}$ and $\phi_{\rm out}$ while avoiding overfitting due to the additional degrees of freedom. One may also adopt a progressive increase of the NA ¹² to utilize the deeper penetration of the forward scattered photons ^{21, 22}, and/or a progressive increase of the depth ⁵⁸. Global optimization can be used when the scattering is so strong that a local optimization loses its consistency. Future work may also explore other optimization schemes beyond using an image quality metric.

Validation against the ground truth is nontrivial for volumetric imaging ⁷⁵. Using the recent “augmented partial factorization” simulation method ⁷⁶, we have validated SMT through full-wave numerical simulations ⁷⁷.

One can interpret RCM, OCT, and OCM as measuring the diagonal elements of the reflection matrix in a spatial basis. SMT utilizes the additional off-diagonal elements to digitally refocus to different planes and to enable double-path wavefront corrections. The price is having to measure those off-diagonal elements.

Motion artifacts make wavefront shaping in vivo a challenge ⁷⁸. While the virtual wavefront shaping of SMT dispenses with slow SLMs, a coherent synthesis still requires measuring the matrix elements before the scatterer arrangement changes substantially. The speckle decorrelation time, primarily limited by the blood flow, is estimated as 5 ms at $\lambda=633$ nm wavelength 1 mm deep inside the mouse brain where blood is rich ⁷⁹ but can go beyond 30 seconds for the skull where there is less blood ⁵⁵. Currently, our total measurement time for the USAF-target-under-tissue sample is three minutes, with the majority of which spent on the mechanical acceleration and deceleration of a birefringent filter during a scan-and-stop operation of the tunable laser. Scanning the filter continuously can reduce the measurement time to 5 seconds (Supplementary Sec. I D). One can further accelerate by orders of magnitude through truncating the spatial measurement ⁵⁷, parallelizing the spectral acquisition ⁸⁰, and reducing the number of frequencies or incident angles (Supplementary Sect. V–VI).

While SMT uses contrasts in the scattering strength, digital staining ^{81, 82} may be used to infer other contrasts. The phase information in $\psi_{\rm SMT}$ can help detect sub-nm displacements ^{83, 84}. One may incorporate polarization gating to select birefringent objects such as directionally oriented tissues. The hyperspectral scattering matrix can additionally resolve spectral information of the sample, such as the oxygenation of the hemoglobin. Applications that require feedback to identify the region of interest can use a subset of the matrix data to reconstruct a coarse image in real time.