1. INTRODUCTION

Due to scattering, high resolution optical imaging of biological tissues is mostly limited to superficial layers. The resolution of fluorescence imaging depends on the ability to form a sharp focus at a plane of interest, which becomes increasingly challenging for tissue layers thicker than a few scattering mean free paths (MFP) [1]. With the advent of wavefront shaping, it became possible to control light propagation and focus both through [2,3] and inside [4–6] turbid media to a diffraction-limited spot. However, this focus spot must be scanned around within the scattering medium to form an image. Assuming the plane of interest is at a depth where excitation light arrives after multiple scattering events, different uncorrelated incident wavefronts are needed to focus the light into neighboring points. This would require a new time-consuming round of wavefront correction for every scanned point, making imaging practically impossible for most applications.

Fortunately, first-order spatial correlations between the incident and scattered wavefronts provide a spatial range where single-wavefront correction remains accurate. This range is referred to as the memory effect [7], with several related varieties [8]. The angular or “tilt/tilt” memory effect describes the following phenomenon: small tilts of the incident wavefront will result in corresponding tilts of the scattered wavefront at the output surface of a scattering medium. This tilt will then manifest itself as a shift of the resulting speckle pattern at a distant imaging plane. Although the “tilt/tilt” memory effect has been exploited in multiple proof-of-principle experiments [9–12], the requirement that the imaging plane is at a distance away from the scatterer renders the tilt/tilt effect of less interest for biological imaging, where scanning must be achieved within tissue at the desired imaging plane. Recently, the translational memory effect—the “shift/shift” memory effect—was discovered and exploited for focusing through forward-scattering media up to one transport mean free path in depth [13]. In this case, a shift of the incident wavefront in the plane orthogonal to the propagation direction results in a shift of the scattered light speckle pattern at the output plane. The “shift/shift” memory effect range was shown to be equal to the average speckle grain size under a plane wave illumination, which provides a convenient way to measure the correlation range, or field of view (FOV). An optimized incident wavefront can then be shifted across the correlation range to scan the focus along the output plane by the same amount.

In the general case of multiple-scattering media, the “shift/shift” memory effect is of limited use for scanning, since under a plane wave illumination light would scatter in a wide output cone, resulting in a small speckle grain size [14]. In biological tissues, however, light scattering is highly anisotropic in the forward direction at visible and near-infrared wavelengths [15]. This means that most photons traveling through tissue would scatter only in a narrow forward cone (so called “snake” photons [16,17]), resulting in long transport mean free paths (e.g., 100 μm MFP and 1 mm TMFP [1]). As snake photons would deviate less from the ballistic trajectory and exit the scattering medium at low angles, we hypothesize that they would exhibit a higher level of “shift/shift” correlations.

Because the time that scattered light spends inside the tissue depends on the scattering path length, weakly scattered snake photons can be separated from multiple-scattered light by means of coherence gating. Coherence gating can be done either in the time [18] or the frequency domain [17], with the frequency domain being preferable, since the spectrum of the time gate can be digitally shaped.

In this work we studied the dependence of speckle spatial correlations on the time light spends in the scattering medium for forward-scattering samples. We reconstructed the time-varying electric field at the output surface of a scattering medium by recording 801 interferograms of speckle patterns with wavelengths ranging from 690 to 940 nm. Performing a Fourier transform along the frequency axis of this dataset, we obtained what is effectively a temporally resolved electric field response to a broadband illumination, with less than 9 fs sampling resolution. We then quantified the extent of speckle translational correlations at different time delays and showed that coherence gating of early arriving photons increases the range of the translational memory effect almost fourfold.

2. METHODS

A. Translational Memory Effect of the Scattered Light Gated in the Frequency Domain

In this section, we will explain how one can form a scattered electric field response to a bandwidth-limited broadband pulse from an array of complex electric fields under monochromatic illumination.

Controlling scattered ultrashort pulses and broadband light in the spectral domain allows for precise manipulation of incident light in a range of applications, including scanningless optical sectioning [19]. A spectral domain approach to time gating was first proposed in Ref. [20,21] and recently used to temporally refocus femtosecond pulses through a scattering medium using the spatio-spectral transmission matrix [22]. A laser pulse can be represented as the sum of its monochromatic components through the Fourier transform. Then, if we record a discrete set of speckle patterns at $N$ frequencies ${\nu }_1\dots {\nu }_N$ with an equal step $\mathrm{\Delta }\nu$, we can approximate the output wavefront as

The set of scattered monochromatic wavefronts at different frequencies ${\{u_{\mathrm{out}}(x,y,{\nu }_n)\}}_{n=1\dots N}$ can be used to reconstruct the scattered wavefronts at different time delays post hoc by performing a Fourier transform along the frequency axis of the dataset. Because of discretization, the frequency step $\mathrm{\Delta }\nu$ together with the spanned frequency bandwidth will define the time periodicity and the temporal resolution of the synthesized pulse [21].

Due to factors such as dispersion, chromatic focus shift, and reference beam phase drift, the recorded monochromatic electric fields are shifted in phase by different unknown global offsets ${\{u_{\mathrm{out}}(x,y,{\nu }_n)e^{i\varphi ({\nu }_n)}\}}_{n=1\dots N}$, where ${\{\varphi ({\nu }_n)\}}_{n=1\dots N}$ are unknown spectral phases [23]. When unaccounted for, these global phase shifts essentially stretch the synthesized pulse in time and degrade the temporal resolution. To try to determine them, we start with the assumption that for a shift in frequency much smaller than the spectral bandwidth of the medium [2], change in the speckle pattern can be considered negligible. We can discretize the output plane into $M\times M$ pixels with indices $k,l=1\dots M$. Reconstructed output fields at neighboring frequencies can then be expressed as $u_{kl}({\nu }_1)=|u_{kl}({\nu }_1)|e^{i{\varphi }_{kl}({\nu }_1)}$ and $u_{kl}({\nu }_2)=|u_{kl}({\nu }_2)|e^{i{\varphi }_{kl}({\nu }_2)}$, each $M\times M$ pixels in size, and with the total intensity normalized to one. If they differ only in global phase shift $\mathrm{\Delta }\varphi$, we can estimate this phase shift by taking an inner product across all pixels,

The argument of the inner product on the left side of Eq. (2) gives us the global phase shift between the two reconstructed electric fields. If we compute and then subtract this phase from the dataset in a chain-like manner (between first and second, second and third frequency measurement etc.), we can align all of the measured monochromatic electric fields at different frequencies to the same global phase.

We can now summarize our method of measuring time-gated scattered fields in the spectral domain [Fig. 2(a)]. We first recorded speckle patterns at a set of equally spaced frequencies, then reconstructed a 3D dataset of scattered electric fields at different frequencies $u_{\mathrm{out}}(x,y,\nu )$, aligned all of the fields to the same global phase via the procedure outlined above, and finally computed the 1D fast Fourier transform (FFT) along the frequency dimension of this 3D dataset. The resulting dataset contains the scattered electric fields at the output plane at different time delays, $u_{\mathrm{out}}(x,y,t)$.

We note that imprecise alignment and accumulation of spectral phase compensation errors led to the remaining linear phase ramp in the dataset, which caused our reconstructed dataset in the time domain to include a global temporal offset. We manually compensated this offset for each scattering sample in order to realign all output pulses to the same initial time delay, which we define as the arrival of ballistic photons at the output surface.

B. Experimental Setup and Sample Preparation

We used a Mach–Zehnder interferometer to record and extract scattered complex electric fields at the output surface of the inhomogeneous medium (Fig. 1). The beam from a Ti:Sapphire laser set to tunable continuous-wave (CW) mode (MaiTai Spectra-Physics) was expanded and attenuated using a half-wave plate and a polarizing beamsplitter (not shown in the figure) before the interferometric setup, resulting in a horizontal polarization of the laser light that enters the interferometric setup. The first beamsplitter (BS1, non-polarizing) separates the light into two arms: the reference arm and the object arm that illuminates the sample. The combination of the microscope objective OBJ2 (Nikon 40x, NA 0.75) and the tube lens L2 ($f=200\mathrm{mm}$) images the scattered light onto an imaging sensor (Basler CMOS2), where it interferes with the reference beam in an off-axis configuration. After switching the laser source to CW mode, we varied the wavelength output from 690 to 940 nm with 801 equally spaced frequency steps and recorded a set of interferograms.

 

Fig. 1. Principal optical setup. The output of a tunable laser source is split into two arms of a Mach–Zehnder interferometer by the beamsplitter BS1. The combination of lens L1 and microscope objective OBJ1 projects a plane wave at the input plane of the scattering medium. A focus spot can be formed and precisely positioned at the same plane by removing L1 and focusing using a reflection imaging system (lens L3 and CMOS1). The scattered light at the output plane of the medium is imaged by the transmission imaging system (OBJ2 and L2) onto the CMOS2, where it is combined with the reference beam, forming a digital hologram. Global phase delay introduced by imaging elements in the object arm is corrected as described in the Methods section. Inset: The translational memory effect can be measured in two ways: (a) as the width of speckle autocorrelation resulting from plane wave input, and (b) as the cross-correlation width between fields produced by physically translating the sample under focused light excitation.

Download Full Size | PDF

Scattered electric fields were reconstructed for each frequency and normalized to unitary power. Weighting different frequencies allows us to shape the spectrum of the illumination source together with its temporal profile. We chose a uniform spectrum shape resulting in a sinc temporal profile with a time-sampling period of $\sim 9\mathrm{fs}$ and a time period of $\sim 7\mathrm{ps}$. According to Eq. (1), a Fourier transform along the frequency axis produced an electric field at different time delays.

We measured the extent of the translational memory effect for different output wavefront delay times in two independent ways (see inset on Fig. 1). As reported in Ref. [13], for an anisotropically scattering medium the translational memory effect range can be measured directly by shifting the input beam across the input surface of the scattering medium and measuring the width of the decorrelation curve for the output speckle pattern. This decorrelation range has been shown to be equivalent to the mean speckle grain size (speckle pattern autocorrelation) that results from illuminating a homogeneous scattering sample with a plane wave.

First, we illuminated the sample with a plane wave projected onto the input surface of the sample by the combination of tube lens L1 and microscope objective OBJ1 (Nikon 20x, NA 0.75), and reconstructed the electric field at different time delays as described above. The range of the translational memory effect was then calculated from the mean speckle grain size at the output surface of the scattering medium, defined as the full width at half-maximum (FWHM) of the speckle autocorrelation. In a second experiment, we computed the translational memory effect range by focusing light onto the input surface of the sample and then physically translating the sample while measuring the cross-correlation of the output electric fields.

Most biological tissues exhibit highly anisotropic forward scattering with a scattering anisotropy factor $g>0.9$ [15]. To mimic the anisotropic behavior of biological tissue, we prepared samples of varying thickness consisting of 5 μm silica microspheres (Sigma-Aldrich) in 1% agarose. Using a Mie calculator [24], the scattering properties of the sample were estimated to be $g=0.976$ for the anisotropy factor and 90 μm for the scattering MFP.

3. RESULTS AND DISCUSSION

Starting with a 360 μm thick sample (equal to four scattering MFP), we reconstructed the temporal evolution of the speckle patterns at the output surface of the scattering medium under plane wave illumination [Fig. 2(b)]. The speckle pattern changed as a function of time, with the size of the speckle grains decreasing with increasing time delay. The average speckle grain size across different time positions was approximately equal to the average speckle grain size produced under CW illumination. After performing a 2D spatial Fourier transform of the resulting complex speckle patterns at different time delays, we observed a ring-like spread of $k$-vectors, with the radius of the ring growing as a function of time [Fig. 2(c)]. This indicates that later-arriving scattered photons emerged from the medium at increasingly larger angles. The experiment was repeated for two thicker samples of 720 and 1080 μm (8 MFP and 12 MFP, respectively). As sample thickness grows, so does the thickness of the rings, since more scattering trajectories can contribute to an output field at a given time delay (Supplement 1, Fig. S1).

 

Fig. 2. Digital time gating of the scattered light in the spectral domain. (a) Sketch of digitally time-gated electric field reconstruction. The amplitudes and phases of the complex speckle fields at the output of a scattering medium are measured for a set of equally spaced frequencies (left), which are then brought to the same global phase (middle). The scatterer’s temporal response is reconstructed by taking a Fourier transform across the frequency dimension (right). (b) Amplitude and phase (color) of the complex speckle pattern at different time delays for the thinnest sample (scale bar, 10 μm). (c) Spatial Fourier transform of each resulting complex speckle pattern (scale bar, absolute value of $\mathit{k}$ that corresponds to a mean single scattering angle, $\sim 13$°). At larger time delays, photons emerge from the medium at an increasing angle. All images were normalized to unitary power.

Download Full Size | PDF

The results of the speckle autocorrelation at different time delays are summarized in Fig. 3(a). The mean speckle grain size was the largest at the earliest time delay, and it quickly decreased at later time delays [Fig. 3(a)]. For three different samples we compared the speckle autocorrelation of the early-arriving scattered light with that of CW illumination at a central wavelength of 815 nm [Fig. 3(b)]. In each sample we observed that the earlier-arriving photons exhibit a wider range of translational memory effect compared to the CW case, as manifested by a larger mean speckle grain size for plane wave illumination. For the thinnest sample, the earliest time delay gave almost a fourfold increase in the memory effect range, while for thicker samples, where ballistic light is virtually absent, it was still more than twofold. Since the synthesis of the pulse was done digitally in the frequency domain, we could select a different bandwidth of the excitation source and therefore create different time-sampling periods. Comparing the results of the memory effect extent as a function of excitation bandwidth [Fig. 4(a)], we observed that the mean speckle grain size (the range of the translational memory effect) increased with larger bandwidths (i.e., with shorter time-sampling period).

 

Fig. 3. Speckle autocorrelation of digitally time-gated light. (a) Speckle autocorrelation as a function of time delay for three samples (under plane wave illumination) with varying thickness, 360 μm (left), 720 μm (center), and 1080 μm (right). (b) Speckle autocorrelation for the earliest time delay compared to the non-gated speckle (CW illumination, 815 nm). The translational memory effect range is wider for earlier-arriving light as compared to CW illumination.

Download Full Size | PDF

 

Fig. 4. Influence of bandwidth on the translational memory effect range and comparison of this range for early-arriving “snake” photons versus CW response based on two independent measurements. (a) Dependence of the speckle autocorrelation FWHM (mean speckle grain size) on the bandwidth of the synthesized pulse for the thinnest sample. A short (broadband) pulse is essential for observing an increased memory effect for early-arriving snake photons. (b) Speckle autocorrelation value (autocorrelation function, ACF) measured for an input plane wave illumination (solid lines) and the “shift/shift” memory effect (cross-correlation function, XCF) correlation value measured for an input focus (dotted lines). Here we show ACF and XCF measurements for both CW illumination (central wavelength 815 nm, blue curves) and for digitally gated early-arriving “snake” photons (red curves). Both measurements indicate an increase in the FOV of the translational memory effect. Data for the early-arriving “snake” photons presented as mean ± standard error of the mean.

Download Full Size | PDF

For the second measurement, we formed a focus spot on the input plane of the scattering sample by removing the tube lens L1, and again recorded a set of electric fields at different frequencies. The focus was positioned at the input plane of the scattering medium with the help of a reflection imaging system (tube lens L3 and Basler CMOS1) placed before the focusing objective OBJ1. We then physically translated the sample in micrometer steps, perpendicular to the incident beam direction and relative to the input focus spot, and repeated the measurement for a total translation distance of 4 μm. The shift/shift correlation function $C_{\mathrm{shift}/\mathrm{shift}}(n\mathrm{\Delta }x)$ was obtained as the absolute value of the cross-correlation between the complex scattered electric fields at two positions with a relative distance of $n\mathrm{\Delta }x$.

Since the shift/shift correlation function only depends on the relative distance between any two points and not their absolute positions, we first created a $m\times m$ matrix of pairwise correlation coefficients $\Vert C_{ij}\Vert$, where $C_{ij}$ is the absolute value of the correlation coefficient between the electric fields at shift positions $i\mathrm{\Delta }x$ and $j\mathrm{\Delta }x$. The final shift/shift correlation function value at a shift can then be calculated as the average along the corresponding side diagonal of the pairwise correlation matrix,

Calculating the correlation coefficient with Eq. (3) allows for more pairwise correlation coefficients to be averaged together for shorter shifts (position closer to the main diagonal of the matrix), therefore offering a better approximation to the shape of the shift/shift correlation function, especially at its onset.

We compared two independent measurements of the translational memory effect range both for the earliest arriving time-gated light and the response under CW excitation (center wavelength of 815 nm, sample thickness 360 μm, MFP 70 μm), in Fig. 4(b). Our two independent measurements of the memory effect range, via speckle autocorrelation and via sample shifting, both matched and confirmed that earlier-arriving light exhibits a larger range of translational correlations.

4. CONCLUSIONS

Correlations in scattered light in the form of the optical memory effect [7] have received considerable attention in the last decades. In its classical formulation, the memory effect requires scattering samples to be thin, which would make it inapplicable to imaging through biological tissues and other thick scattering media. However, recent studies have shown that even thick ($\sim 1\mathrm{mm}$) biological media exhibit strong scattering correlations [13,25]. Newly discovered “shift/shift” correlations extend deep inside the scattering tissue, as long as scattering is anisotropic and preserves some degree of directionality between the input and output wavefronts. In theory, this effect should be the strongest for snake photons that are deflected at small angles from their original direction. However, under CW illumination those photons cannot be separated from multiply scattered light, so their properties cannot be exploited.

In this study we digitally synthesized an ultra-short pulse from a set of monochromatic holograms recorded at different frequencies. Broadband illumination allowed us to reliably gate only the early-arriving photons and demonstrate that the range of the translational memory effect that those photons exhibit is almost fourfold larger than the average effect under CW illumination. This FOV increase was clearly observed for all studied scattering samples, with the FOV decreasing slowly as the sample thickness increased and the proportion of weakly scattered light diminished.

We employed most of the laser bandwidth that could be addressed for uninterrupted wavelength tuning. As demonstrated by the dependency of the translational memory effect range on the illumination bandwidth [Fig. 4(a)], a wide bandwidth was essential for achieving high enough temporal resolution to reliably select only early-arriving light. This can be explained by the fact that highly anisotropic forward-scattering media have short confinement times, and the emergent output pulse is relatively weakly elongated compared to the incident pulse.

Time-gating of early-arriving “snake” photons with wide bandwidth sources could be combined with nonlinear deep tissue imaging methods that rely on wavefront correction to focus through scattering tissue layers [26–28]. These techniques indirectly select the first-arriving light through iterative wavefront optimization. Adding the temporal degree of freedom in wavefront shaping can enable more precise control of light propagation through scattering media, therefore paving the way towards more efficient scattering compensation for deep tissue imaging.