Optical imaging plays a crucial role in nondestructive, high-resolution visualization, owing to its high spatial resolution. It has been widely employed to study biological and physical processes, enabling applications such as cellular-level neuronal activity monitoring, diagnostic imaging in ophthalmology and audiology, and surgery-free biopsies. Beyond biology, optical imaging also holds great promise for applications such as optical characterization of bulk colloids and high-resolution, nondestructive device testing. However, optical imaging is significantly limited by its inability to penetrate deeply into complex scattering media, such as biological tissue, where inhomogeneities lead to aberrations and multiple scattering that degrade image quality even to the point of being unrecognizable.1 Although aberrations distort the light wavefront and reduce spatial resolving power, severely undermining the image contrast, multiple scattering exponentially attenuates the signal and introduces a speckled background, obscuring the already-aberrated desired information.2^,3 Overcoming these challenges to achieve label-free, high-resolution imaging through complex scattering media is essential to expand the utility of optical imaging in practical applications.

Over the years, gating techniques have been developed to enhance image contrast, resolution, and depth in optical imaging by selectively isolating the single-scattered light from the unwanted background. For example, reflectance confocal microscopy (RCM)4^,5 uses spatial gating. Optical coherence tomography (OCT)6 and optical coherence microscopy (OCM)7 apply a coherence time gate in addition. Other approaches include selective detection of the forward scattered photons8^,9 and photoacoustic tomography,10 which image deeper with a reduced resolution. Despite these efforts, these techniques are generally limited to a depth-over-resolution ratio of around or below 200 in biological tissue10^–12 due to the inherent exponential decay of single-scattered light, meaning the high resolution is confined to superficial layers. To further complicate matters, the single-scattered wave is affected by the unavoidable aberrations induced by the sample’s heterogeneity, severely degrading image quality.

Adaptive optics (AO) has been a significant advancement in improving the quality of the spatial gate by correcting for low-order wavefront distortions13^,14 within a volume called the isoplanatic patch.14^,15 Traditional AO systems rely on guide stars or wavefront sensors to correct aberrations, whereas the image-based sensorless AO16 using Zernike-mode-based image metric optimization eliminates this need but still uses hardware with limited correction modes and speed. Computational AO17^–20 offers potentially more correction modes, faster speed, and better multiple-scattering immunity. Despite these advancements, AO falls short of achieving high-resolution imaging in the presence of strong high-order aberrations and scattering. Higher-order aberrations require correcting both the incident and the outgoing wavefronts,21 but existing computational AO schemes can only perform one correction. Furthermore, the severe multiple scattering causes large errors in AO, leading to a breakdown in performance and limiting its applicability in deeply scattering media.

Recently, the matrix-based approach22^,23 has gained significant attention for its potential to push the depth limits of optical imaging using the full sample’s optical response. Unlike confocal detection, the matrix-based approach captures signals arriving at both confocal and nonconfocal positions, building a more complete picture of the sample’s input–output responses, known as the scattering matrix,24^,25 which allows for the digital correction of both the incident and reflected wavefronts. Pioneering works on the “closed-loop accumulation of single-scattering” (CLASS) method26^–30 and the “distortion matrix” approach31^–33 (whose latest version is equivalent to CLASS34) have leveraged correlations within the reflection matrix, the reflection mode of the scattering matrix, to correct wavefront distortions. Extending to the multispectral matrix imaging scheme,33^,35 the recent “volumetric reflection-matrix microscopy” (VRM)35 generalizes CLASS to add dispersion compensation. However, current matrix-based methods still face significant challenges when the single scattering signal is severely distorted by high-order sample-induced aberrations and overwhelmed by the multiple scattering backgrounds, which introduce significant errors in the wavefront correction. To fully exploit the wealth of information contained in the scattering matrix, robust and effective wavefront correction approaches are required—ones that can withstand the detrimental effect of the strong multiple scattering to undo aberrations and unlock more potential of matrix-based imaging for high-resolution imaging through complex media.

Here, we propose an adaptive optical multispectral matrix approach to tackle the challenge of achieving high-resolution, label-free tomographic imaging through complex scattering media. Relying on the sample’s full far-field linear response encapsulated in the scattering matrix, the approach, termed scattering matrix tomography (SMT), serves as a unified platform for systematically exploring the rich information embedded in the scattering matrix, allowing us to perform scattering-based imaging with various gating mechanisms and apply flexible corrections. Furthermore, inspired by AO, we formulate the problem of imaging through complex scattering media as a numerical optimization rather than relying on matrix correlations. To overcome the significant challenges posed by high dimensionality and nonconvexity due to severe sample-induced aberrations and optimization error due to overwhelming multiple scattering, we employ wavefront regularizations and coarse-to-fine progressive optimization strategy along with additional refinements that further enhance the robustness and image quality. Together, these advancements enable our method to achieve 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 nonoptical label-free methods) to our knowledge—where the signal is reduced by over 10-million-fold due to multiple scattering and is completely overwhelmed by speckles prior to the digital corrections. 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. We synthesize conventional reflection-based methods (RCM, OCT, and OCM) and the most advanced matrix-based method VRM and found all of them to fail at these depths. By framing the challenges of imaging through scattering media as numerical reconstruction and optimization problems, employing wavefront regularization and the progressive optimization strategy, our method builds upon the strengths of AO and the scattering matrix while providing important advancements that enable deeper, higher-resolution imaging in complex scattering environments, offering an intuitive and versatile platform that uses wave physics and computation to push the capabilities of imaging.

The scattering matrix S(kout,kin,ω) encapsulates the sample’s far-field linear response24^,25 [Fig. 1(a)]. At frequency ω, for each incidence with momentum kin, the scattered fields are measured and stored at a corresponding matrix column, with each row corresponding to one output momentum kout. As any wavefront can be synthesized from plane waves with designed phases and amplitudes, the scattering matrix, with its flexible wavefront engineering capability, allows one to mimic imaging experiments with any arbitrary incident and scattered waves by computationally tailoring the matrix columns and rows, thus achieving desirable images.

In light of this, as illustrated in Figs. 1(b) and 1(c), after a one-time measurement of a subset of S(kout,kin,ω), we utilize this data to perform virtual imaging experiments—digitally synthesizing the would-be response of the sample for any customized spatiotemporal input and any tailored measurement in space-time without additional physical experiments. In these virtual experiments, we can enforce all gating mechanisms (time gating and confocal spatial gating) and all possible correction schemes (pulsing shaping, index-mismatch correction, and double-path corrections of high-order aberrations), scan the focus, and optimize the image quality, all carried out digitally without having to wait for the hardware. Below, we provide a rigorous framework to reconstruct images and correct for wavefront distortions using the scattering matrix.

Any incident wave is a superposition of plane waves: Ein(r,t)=∑kin,ωE˜in(kin,ω)eikin·r−iωt, where r=(x,y,z) is the position, t is time, and E˜in(kin,ω) is the amplitude of its plane-wave component with incoming momentum kin at frequency ω. The amplitudes E˜out(kout,ω) that form the resulting outgoing wave Eout(r,t)=∑kout,ωE˜out(kout,ω)eikout·r−iωt are given by the scattering matrix through E˜out(kout,ω)=∑kinS(kout,kin,ω)E˜in(kin,ω). The angular summations are restricted to |k|2=|k∥|2+kz2=(nbgω/c)2 in a background medium with speed of light c/nbg; we use “angle” interchangeably with “momentum.”

We perform the digital gatings through summations [Fig. 1(d)]: summing over plane waves eikin·r with incident angle kin and amplitude E˜in(kin,ω)=e−ikin·rin creates an incident beam Ein(r,ω) spatially focused at rin, forming an input spatial gate; summing over plane waves eikout·rout with outgoing angle kout and amplitude E˜out(kout,ω) from the scattering matrix yields the scattered field Eout(r=rout,ω) given a virtual detector at a position conjugated to rout, forming an output spatial gate; summing over frequencies ω with an e−iω(t−t0) phase creates an incident pulse that arrives at rin at time t=t0, forming a temporal gate. Given a spatiotemporal focus arriving at rin at time t0, the scattered field at position rout at time t is therefore Eout(rout,t)=∑kout,kin,ωS(kout,kin,ω)eikout·rout−ikin·rin−iω(t−t0). In our virtual experiment, we align the input spatial focus with the output spatial detection (setting rin=rout=r) and align the temporal focus with the confocal spatial one (evaluating Eout at time t=t0), which yields the triply-gated scattering amplitude of the sample at position r, denoted as ψSMT(r)=∑kout,kin,ωei(kout−kin)·rS(kout,kin,ω).(1)

Digitally scanning r [Fig. 1(c)] 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 kout−kin in the summation (typically set by the numerical aperture 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. Like digital holography,36 the ψSMT(r) image can be volumetric or a slice with any orientation. 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.37 One can utilize any subset of the scattering matrix—reflection, transmission, remission,8^,38 a combination of them, or other subsets—with any number of angles and frequencies. We use the nonuniform fast Fourier transform39 to efficiently evaluate these summations and the 3D spatial scan. Equation (1) is the minimal version of SMT.

The triple gating suppresses multiple scattering. The single-scattering field at output kout for an input from kin is given by the Born approximation as proportional to η˜(q), namely, the q=kout−kin component of the 3D Fourier transform of the sample’s permittivity contrast profile η(r).6^,40^–42 Such single-scattering contributions add up in phase in Eq. (1) to form the image, similar to an inverse Fourier transform from η˜(q) to η(r). The multiple-scattering contributions do not add up in phase as they differ strongly among different kout, kin, and ω. In particular, coherently summing over N inputs, outputs, and frequencies raises the single scattering signals N times, whereas the multiple scattering signals only go up by N times due to their incoherence.22 Therefore, the triple summations over kout, kin, and ω boost the single-to-multiple-scattering ratio in ψSMT(r). From this reasoning, in theory, as long as there are enough inputs, outputs, and frequencies, multiple scattering can be suppressed below single scattering, leading to deeper imaging depth.

Next, we perform spatiotemporal corrections in the virtual experiments to reverse aberrations by computationally engineering the phases of the incident and scattered wavefronts. 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. We compensate for such dispersion using a spectral phase θ(ω) acting as a virtual pulse shaper [Fig. 1(e)], such as in spectral-domain OCT.18^,43 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,44 and misaligns the spatial gates and the temporal gate. Using the appropriate propagation phase shift (kout−kin)·r in each medium and the appropriate arrival time, we restore 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. 1(f)] without liquid immersion (Sec. III C in the Supplementary Material). Importantly, we also digitally introduce angle-dependent phase profiles ϕin(kin) and ϕout(kout) that act as two virtual spatial light modulators (SLMs), one for the incident wave and one for the outgoing wave [Fig. 1(g)], to compensate for aberrations. These two virtual SLMs can take complicated phase profiles and thus are able to perform extreme cases of AO where the corrected sample-induced aberrations are of high order and fast-varying. This yields SMT with corrections $$\begin{gathered} ISMT(r)=|∑ωeiθ(ω)∑kouteikout·r+iϕout(kout) \\ ⁢∑kine−ikin·r+iϕin(kin)S(kout,kin,ω)|2. \end{gathered}$$(2)By measuring the reflection from a mirror, we remove the dispersion and aberrations in the input path of the optical system (Sec. III B in the Supplementary Material).

Next, we look for suitable corrections θ(ω), ϕin(kin), and ϕout(kout). Enhancing the focus at r will enhance the scattering intensity from a target there, so ISMT(r) acts as an intrinsic virtual guidestar at r [Fig. 1(c)]. Therefore, we maximize an image quality metric M, {θ,ϕin,ϕout}=argmaxθ,ϕin,ϕoutM,M=∑rISMT(r)lnISMT(r)I0,(3)with ∑r being a volumetric summation and I0 being a constant; the logarithm factor promotes image sharpness.45 Similar image metrics have been used in single-path computational AO18 and single-path isoplanatic wavefront shaping with an SLM46^,47; here, we adopt it for digital double-path spatiotemporal corrections. We derive the gradient of M with respect to the parameters and use a quasi-Newton algorithm, the low-storage Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method,48 in the NLopt library49 to maximize M. Since a wavefront correction ϕin(kin) or ϕout(kout) only remains optimal within an isoplanatic patch,14^,50^,51 we divide the space into zones (both axially and laterally) and use different ϕ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 [Figs. 2(a)–2(f)]. This happens because (1) the problem is high-dimensional and nonconvex with numerous poor local optima, (2) the large number of parameters in θ(ω), ϕin(kin), and ϕout(kout) in the many zones can lead to overfitting, where the optimization inadvertently enhances unintended contributions like the speckled background instead of signals from a target at the intended position r, and (3) the strong quasi-random multiple scattering signals lead to phase errors in the wavefront correction. To mitigate these issues, we employ five strategies to guide the search, avoiding the poor local optima, overfitting, and the phase errors caused by multiple scattering (Sec. III E and F in the Supplementary Material). (i) Regularize θ(ω) to a third-order polynomial43 to avoid overfitting. (ii) Regularize ϕin(kin) and ϕout(kout) each to a Zernike polynomial14^,16 with a controlled number of terms to avoid overfitting and errors by multiple scattering. (iii) First optimize ϕin/out over the full volume and then progressively optimize over smaller zones while using the previous ϕin/out as the initial guess to efficiently exploit the spatial correlation of sample-induced aberrations, avoiding suboptimal local optima. (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 because the low-order terms correct for slowly varying aberrations with a large isoplanatic volume, whereas the high-order terms correct for the fast varying aberrations with a small isoplanatic volume. With these strategies, we find corrections θ(ω), ϕin(kin), and ϕout(kout) that clearly reveal the target [Figs. 2(g)–2(k)] even though the initial image was overwhelmed by speckles due to strong multiple scattering despite triple gating [Fig. 2(a)].

Although triple gating by itself is not sufficient, it is a necessary ingredient in the optimization because it promotes the target-to-speckle ratio. Without the summations over kout, kin, and ω in Eq. (2), the optimization may maximize the background speckles instead when scattering is strong. As all targets within an isoplanatic volume share the same optimal wavefront, whereas the speckled backgrounds do not, the volumetric summation over space r in Eq. (3) also promotes the target-to-speckle weight in the gradient of the optimization. Figure S10 in the Supplementary Material shows that without one of these quadruple summations in {r,kout,kin,ω}, the optimization landscape may not reflect the target contribution and may exhibit many local optima associated with the multiple-scattering background. This suggests that optimizing over a spatial region where the wavefront distortion remains approximately constant, that is, within an isoplanatic area, is important to improve the optimization robustness and that the optimization will be less affected by poor local optima if the isoplanatic patch is large. A notable example is the recently proposed image-guided computational holographic wavefront shaping method by Haim et al.,52 which successfully demonstrated detection-path correction across a large number of scattered modes using iterative wavefront optimization. Although this work achieved high-quality imaging of targets located at a distance behind the complex medium, many real-world applications involve more challenging conditions, such as targets embedded within volumetric scattering media where both the illumination and detection paths are affected by strong, spatially varying distortions. These conditions result in much smaller isoplanatic volumes and stronger multiple scattering, making direct optimization prone to overfitting or convergence to poor local optima and speckle-dominated solutions. SMT addresses such challenges by conducting double-path corrections along with regularization and coarse-to-fine optimization strategies.

CLASS26^,28^–30 (which is equivalent to the distortion matrix approach32^–34) and VRM35 are based on the correlations of the single-scattering signal; they were not formulated as a maximization problem like Eq. (3). However, our analysis reveals that they in fact implicitly maximize some metrics (Sec. IV D in the Supplementary Material). However, these implicit maximizations were not designed to avoid poor local optima or overfitting, which are important (Fig. 2). In addition, the dispersion compensation phase for each wavelength of VRM35 is found by correlating the reflected fields of each input at that wavelength and the central wavelength, which are without a confocal gate, a temporal gate, and a volumetric spatial summation, all of which are important (Fig. S10 in the Supplementary Material) to reduce the error caused by multiple scattering; the other methods, meanwhile, do not compensate for dispersion. These are part of the reasons why the previous matrix-based methods fail when multiple scattering is strong enough to cause errors in the aberration correction.

We use off-axis holography53 to measure the scattering matrix in reflection mode, for which S(kout,kin,ω) reduces to the reflection matrix R(kout,kin,ω) (Fig. 3). 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 kin, and a tunable laser (M Squared SolsTiS 1600) scans the frequency ω. The data acquisition here is 2 to 3 orders of magnitude faster than previous measurements of broadband scattering matrices.25^,35^,54 We use a dry objective lens (Mitutoyo M Plan Apo NIR 100×, NA = 0.5). The power onto the sample ranges from 0.02 to 0.2 mW. We measure around 250 wavelengths λ from 740 to 940 nm, 2900 outgoing angles, and 3900 incident angles within the NA, over a 51μm×51μm area. The detection sensitivity, currently limited by the residual reflection from the objective lens, is 90 dB. See Secs. I and II in the Supplementary Material for details.

We first image the eighth group (with a sixth element that has a bar width and separation of 1.1μm) of a 1951 USAF resolution target underneath a 0.98-mm-thick tissue slice from the cerebral cortex of a mouse brain [Fig. 4(a)]. A standard bright-field microscope image (with incoherent white-light illumination) shows no feature [Fig. 4(b)].

Before considering SMT, we first use the measured hyperspectral reflection matrix to conduct virtual experiments that mimic existing reflection-based imaging methods (Sec. IV in the Supplementary Material): restricting the frequency summation of Eq. (2) to one frequency yields a synthetic RCM image without time gating [Fig. 4(c)]; restricting the angular summations of Eq. (2) to small angles (we use NA = 0.1 here) and fixing the depth of the spatial focus yields a synthetic OCT image [Fig. 4(d)]; for OCM [Fig. 4(e)], we sum over all frequencies and all angles. The depth of the focal plane is chosen to maximize the total signal. All of these images include corrections for the input aberrations of the optical system. For OCT and OCM, we also perform digital dispersion compensation (following the same steps as in SMT). Despite the confocal spatial gate, temporal gate, and corrections, none of these methods reveal any group-8 element due to the overwhelming scattering from the tissue. To compare to the most advanced matrix-based method, we also implement VRM35 (Sec. IV D) in the Supplementary Material, which performs [Fig. 4(f)] worse than SMT without the regularization and progression strategies [Fig. 2(c)] for reasons explained earlier. In addition, we benchmarked against the eigenchannel-based approach23^,55^,56 with the details and results in Sec. IV E in the Supplementary Material. In this approach, the scattering matrix is decomposed into singular values and vectors corresponding to distinct optical channels, and then, the smaller singular values that correspond to multiple scattering eigenchannels are set to zeros to suppress multiple scattering. Although this method offers better computational efficiency through its single-step decomposition process, our comparative analysis reveals its limitations in imaging performance. Although it reduces multiple scattering effects in certain regions, it fails in most of the fields of view due to the lack of an effective wavefront correction.

In comparison, the SMT image resolves the eighth group of the USAF target with near perfection down to the smallest sixth element [Fig. 4(g)]. Here, the progression goes down to 8×8 zones, with up to 275 Zernike polynomials (22 radial orders) per zone in both ϕin and ϕout. The optimized correction patterns ϕin(k∥in) and ϕout(k∥out) are shown in Fig. 4(h). Here, ϕout(k∥)≠ϕin(−k∥) because the former includes the optical-system aberration. Notably, SMT achieves good reconstructions even when the number of frequencies or incident angles is further reduced by over an order of magnitude, including when the image remains fully speckled after dispersion compensation (Sec. V in the Supplementary Material). The total number of wavefront correction variables, 2×275×82=35,200, is large. Figure S13 in the Supplementary Material shows that all of these terms work together to raise the image metric; there is no redundancy. The wavefront corrections vary faster at larger angles due to the use of Zernike polynomials as wavefront regularizers, which causes larger momenta to have larger correction phases and become fast-varying after phase wrapping. Future works will explore various strategies to correct aberrations at large angles more effectively, such as defining the Zernike polynomials with a larger pupil diameter, alleviating this fast variation.

The phase-resolved interferometric data in the synthetic OCT and OCM can be used for dispersion compensation and digital aberration correction via computational AO.18^–21 However, the correction only applies to the phase image itself; it cannot be applied separately to the incident path and the return path because the reflection matrix is not available in OCT/OCM. Figures 5(a)–5(c) show correction on an OCM image without input spatial gate (similar to Ref. 18), which fails here given the lack of triple gating. Figures 5(d)–5(f) show correction in the reciprocal space of a confocal OCM image (mimicking Refs. 19–21; see Sec. IV C in the Supplementary Material for details), which yields negligible enhancements because most of the wavefront distortions here are high-order and can only be corrected in a double-path configuration.21 Unlike computational AO, SMT enables triple gating and double-path wavefront correction, which is necessary to overcome the strong scattering here [Fig. 5(g)].

To quantify the imaging performance, we obtain the point spread function (PSF) by evaluating Eq. (2) with a variable output position r=rin+Δr given a fixed input position rin on the sixth element: PSF(r)=I(r,rin); see Sec. VI in the Supplementary Material. The PSFs of RCM, OCT, OCM, and VRM [Figs. 4(i)–4(l)] have no discernible peak near Δr=0, consistent with their complete failure to image here. The speckled OCM PSF averages to be 70 dB below the peak PSF of a mirror without the brain tissue [Fig. 4(n)], so the signal (which is buried beneath the speckled background and not visible here) has been reduced by at least 10-million-fold due to multiple scattering. The PSF of SMT [Fig. 4(m)] exhibits a sharp peak at Δr=0 with seven times the height of the tallest speckle in the background, showing SMT has not reached its depth limit (where the signal strength equals the background strength) yet. The SMT peak’s full width at half maximum (FWHM) is 1.08μm, close to the mirror PSF’s 0.93μm FWHM, demonstrating diffraction-limited resolution despite the overwhelming multiple scattering. To our knowledge, the depth-over-FWHM ratio of 910 here is the highest reported in the literature for imaging high-contrast targets inside or behind ex vivo biological tissue, which scatters about twice as much as in vivo tissue.57

We next perform 3D tomography of a dense colloid consisting of high-index titanium dioxide (TiO2, refractive index 2.5) nanoparticles dispersed in polydimethylsiloxane (PDMS, refractive index 1.4). The nanoparticles (Sigma-Aldrich, 914320, St. Louis, Missouri, United States) have a typical diameter of 500 nm, and we use Mie theory to estimate the scattering and transport mean free paths to be ℓsca=0.19mm and ℓtr=0.47mm at λ=840nm (Sec. VII in the Supplementary Material). We measure the reflection matrix over a 51μm×51μm field of view with the reference plane at depth z0=1.525mm. Using the reflection matrix, we reconstruct an SMT image over a 50μm×50μm×110μm volume inside the colloid by dividing the depth of field (DOF) Δz=110μm into 16 sub-volumes vertically.

SMT creates a detailed 3D image of all the nanoparticles in this volume [Fig. 6(a)], with zoom-ins and cross sections of individual particles shown in Figs. 6(f), 6(k), 6(l). A statistical analysis of these particles finds the lateral FWHM to vary from 0.7 to 1.0μm within the 110μm DOF and the axial FWHM very close to the theoretical estimate of 1.42μm given the spectral bandwidth of Δλ=187nm here (Sec. VIII in the Supplementary Material). Figures 6(m) and 6(n) show examples of SMT images of nearby particle pairs, with one pair separated horizontally by 0.94μm and one pair separated vertically by 1.48μm. The SMT DOF covers approximately the volume of overlap among the input/output beams in the scattering matrix measurement; it grows with the field of view and is not restricted by the Rayleigh range zR (Sec. VIII B in the Supplementary Material). Here, zR=nsamλ0/(πNA2)≈1.5μm (with center wavelength λ0=840nm and nsam=1.4 for PDMS), whereas the SMT DOF is 73 times larger.

We also construct VRM (generalized from isolated fixed-depth slices35 to volumetric and vertical slices) and synthetic OCM, OCT, and RCM images [Figs. 6(b)–6(e)] for comparison, with longitudinal and transverse slices shown in Figs. 6(g)–6(j), 6(l). For OCM and OCT, the depth of the spatial gate is fixed on one focal plane at z0=1.525mm when the temporal gate is scanned in z (Sec. IV B in the Supplementary Material). All of these methods fail and cannot identify any of the nanoparticles at these depths due to overwhelming multiple scattering. The initial SMT image is also dominated by speckles (Sec. III F.7 in the Supplementary Material) but eventually reveals nanoparticles with high contrast and resolution thanks to corrections. The VRM image exhibits bright artifacts at several constant-z slices because its dispersion compensation maximizes the intensity at those isolated slices (Sec. IV D in the Supplementary Material).

SMT provides an intuitive and versatile framework for noninvasive label-free computational imaging deep inside scattering media. By formulating the wavefront correction and image reconstruction into a nonconvex optimization problem and regularizing the wavefronts, SMT sets a framework for future matrix-based imaging works to build upon, as the optimization algorithm, the optimized image quality metrics, and the number of optimized wavefront modes, can be flexibly chosen depending on each particular sample. Although only the ex vivo mouse brain and dense nanoparticle colloid are used in this paper to primarily demonstrate the capabilities and potentials of SMT, future works will deal with different types of biological tissues such as bone or muscle and nonbiological scattering media like white paint or white papers as we work to implement SMT in various practical applications such as medical imaging or material inspection.

Compared with existing state-of-the-art matrix-based imaging techniques,28^,33^,35 SMT differs in several key aspects: the quadruple summation in image quality metric greatly boosts target-to-speckle ratio; the coarse-to-fine progressive multiscale correction further mitigates local minima caused by strong multiple scattering background noise; and the Zernike-regularized wavefront parameterization avoids overfitting to prevent the optimization from strengthening speckles rather than the target signal. These strategies enable SMT to perform more robust correction in highly scattering, inhomogeneous samples, pushing the imaging depth and quality beyond what current reflection matrix-based techniques can achieve.

One 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 on the frequency.58 Future work can explore optimization strategies that incorporate the frequency dependence in ϕin and ϕout while avoiding overfitting and the larger number of poor local optima due to the additional degrees of freedom (Sec. III F.8 in the Supplementary Material). One may also explore other optimization schemes beyond using an image quality metric.

Validation against the ground truth is nontrivial for volumetric imaging.59 Using the recent “augmented partial factorization” (APF) simulation method,60 we have validated SMT through full-wave numerical simulations.61 In addition, numerical methods such as APF or Monte Carlo simulations62 may help estimate how many inputs, outputs, and frequencies are needed to effectively suppress multiple scattering below single scattering. In particular, simulations of phantoms provide information on how many photons are scattered multiple times versus going ballistic,9 leading to the ratio ρ between single and multiple scattering signals at a particular imaging depth in real samples. Coherently summing over N inputs, outputs, and frequencies increases the ratio between single and multiple scattering signals by N times.22 Guided by simulations, the number of inputs, outputs, and frequencies enough to suppress multiple scattering should be larger than 1/ρ2. Furthermore, the electric field Monte Carlo methods63^–66 and full-wave simulations such as APF67 can model the propagation of complex electric fields inside scattering media. This could offer valuable physical insights to better understand and characterize multiple scattering signals in terms of field-level properties, such as the electric field amplitude, polarization, and phase shifts, showing great potential to further support and inform the developments of SMT. Numerical simulations of modeled phantoms will therefore be employed in future works to guide the design of SMT measurement setups and to develop advanced multiple scattering suppression strategies.

To ensure successful practical applications, particularly for real-time and in vivo imaging, several factors must be addressed, and further improvements can be made in subsequent steps. One can interpret RCM, OCT, and OCM as measuring the diagonal elements of the reflection matrix in a spatial basis, whereas SMT, as a matrix-based imaging approach, utilizes the additional off-diagonal elements to digitally refocus to different depths and to enable double-path wavefront corrections. The price is having to measure those off-diagonal elements. Although the use of all matrix elements enables high-quality imaging, it presents scalability challenges for large field-of-view applications, particularly for highly dynamic samples. This is because the number of matrix elements scales quadratically with the number of diffraction limits within the field of view. These additional measurements require faster acquisition, especially in in vivo imaging, where motion artifacts make aberration correction and imaging in vivo a challenge.68 Although the virtual imaging and wavefront correction of SMT dispense with slow SLMs, a coherent synthesis still requires us to complete the measurement of the matrix elements before the scatterer arrangement changes substantially. In addition, as a computational imaging and aberration correction method, SMT requires accurate phase information. The motion of samples, especially living samples, can create a significant phase instability. To avoid motion artifacts and motion-induced phase instability, the most straightforward way is to expedite the measurement. The speckle decorrelation time, primarily limited by the blood flow, is estimated as 5 ms at λ=633nm wavelength 1 mm deep inside the mouse brain where blood is rich69 but can go beyond 30 s for the skull where there is less blood.29

These challenges in scalability and in vivo measurements open up important, promising directions. Of particular interest is investigating the relative significance of different matrix elements, which could enable selective measurements to reduce both acquisition time and computational demands. Specifically, instead of measuring the scattering matrix on the angular basis, measurement can be done on the spatial basis where only the reflection from output locations rout close to the incident location rin needs to be collected, dramatically reducing the number of matrix elements to measure,32 leading to faster measurement. In addition, encouragingly, other ongoing efforts to improve the measurement speed of matrix-based imaging,32^,33^,35 alongside advancements in hardware development, are making SMT increasingly practical for real-world applications. Currently, our total measurement time for the USAF-target-under-tissue sample is 3 min, the majority of which was spent on the mechanical acceleration and deceleration of a birefringent filter during a scan-and-stop operation of the tunable laser. Scanning the frequency continuously can reduce the measurement time to 5 s (Sec. I D in the Supplementary Material). One can further accelerate by orders of magnitude by reducing the number of frequencies or incident angles (Sec. V in the Supplementary Material). Advancements in hardware are also playing a key role in improving measurement speed. The use of high-speed swept sources can significantly reduce spectral scanning times, whereas innovations to parallelize the spectral acquisition, such as imaging mapping spectrometers,70 allow for the single-shot capture of multiple spectral images, further expediting the data collection. With such accelerations, in vivo SMT data acquisition is possible even inside blood-rich areas such as the brain.

Beyond accelerating data acquisition, reducing the time spent on computational processing is equally important for practical use, especially for real-time applications. For SMT as well as other matrix-based imaging methods, the primary challenges stem from both data transfer speeds and computational intensity. Although the computations are performed after all data is acquired, so in principle, it does not hinder the prospect of in vivo imaging, the long processing time is a significant challenge for real-time imaging, where all the data must be transferred and processed in a short amount of time. The current implementation of SMT employs nonuniform fast Fourier transform for image reconstruction, achieving competitive processing efficiency compared with conventional matrix multiplication methods common in other matrix-based techniques.32^–34 However, significant challenges remain for real-time applications. Currently, the time spent on processing raw data, finding suitable wavefront correction phases, and reconstructing images may range from a few minutes to hours, whereas data transferring between hardware components can take up to 20 min (Sec. III H in the Supplementary Material) due to the large data size. Potential graphics processing unit (GPU) acceleration,33^,34 coupled with ongoing developments in data transmission technology and camera interfaces, as well as selective matrix element measurements, suggests promising pathways for enhancing SMT’s practical utility, particularly in applications where rapid feedback is crucial.

Another aspect of SMT worth exploring is the ability to integrate it with other imaging modalities for practical applications. When measuring the matrix on an angular basis, SMT illuminates the samples with plane waves with varied incident angles and scanned frequencies. This makes SMT a more generalized version of full-field swept-source OCT.33 Alternatively, SMT can also measure the matrix on a spatial basis, where it employs point illuminations such as in OCT/OCM but captures the scattering from not only the illuminated locations but also from other surrounding positions. Therefore, SMT can be modified to be used in conjunction with OCT/OCM, sharing the same illumination part. In addition, as SMT provides a powerful aberration correction, its correction maps can be useful to counteract aberrations in fluorescence or nonlinear microscopy, where the phase information is usually difficult to obtain, making computational aberration correction difficult. In Sec. IX in the Supplementary Material, we propose some future optical systems where SMT can be integrated into optical-fiber-based swept-source OCT for surgical applications and fluorescence microscopy for biological applications.

Although SMT uses the scattering strength as the contrast, digital staining71^,72 and dynamic variation73^,74 may be used to infer other contrasts. Inspired by OCT angiography,74 SMT can also image the decorrelation of highly dynamic samples such as blood vessels (Sec. X in the Supplementary Material), with the difference being that each SMT image has higher quality than each OCT image due to better aberration correction, potentially revealing more sample features. The phase information in ψSMT can help detect sub-nanometer displacements.75^,76 One may incorporate polarization gating to select birefringent objects such as directionally oriented tissues. As SMT is capable of measuring the full vectorial scattering matrices (Sec. XI in the Supplementary Material), future works can explore how samples interact with light of different polarization. The hyperspectral scattering matrix can additionally resolve spectral information of the sample, such as the oxygenation of the hemoglobin.

Yiwen Zhang received her BS and PhD degrees in condensed matter physics from Fudan University in 2014 and 2019, respectively, investigating nanophotonic materials and advanced optical characterization techniques. From 2019 to 2025, she was a postdoctoral fellow in Ming Hsieh Department of Electrical and Computer Engineering at University of Southern California, focusing on optics in complex media, computational imaging, and wave propagation in multi-mode fibers.

Minh Dinh received the BS degree in control and automation engineering from Hanoi University of Science and Technology in 2021. Since 2022, he has been working toward his PhD in electrical engineering at University of Southern California, focusing on optical imaging for biomedical applications.

Biographies of the other authors are not available.