1. Introduction

Optical imaging has provided non-invasive means of interrogating biological subjects, making it a vital tool to elucidate various biological processes and disease progressions in their most native states. However, its working depth has been so shallow that many important reactions occurring deep within living tissues have been out of reach [1]. Heterogeneous spatial refractive index distribution in tissues induces multiple light scattering, which leads to the loss of image information as the propagating waves randomly change their directions on their way to and from the object of interest embedded within the scattering medium. For this reason, most existing high-resolution imaging modalities have relied on detecting ballistic wave (BW) that preserves propagation direction in the scattering medium. However, the amount of BW decreases exponentially with depth, making the image information readily obscured by multiple-scattered wave (MW) even at shallow depth. Therefore, achievable imaging depth for the diffraction-limited optical imaging is mainly determined by the performance of gating operations that selectively detects BW in the background MW.

Finding BW among the strong MW background has been carried out in two imaging geometries: transmission and reflection. In the transmission geometry, researchers have focused on detecting early-arriving photons in a diffusive regime [2–4]. In the scattering medium, BW passes straight through the medium with shortest pathway, thereby arriving earliest at the detector. On the contrary, the optical path length of MW is elongated as its route is deviated from the straight path by the scattering events [5]. In the diffusive regime, where the scattering medium is thicker than approximately 10 times the transport mean free path, BW is well separated with MW in a time domain. Therefore, either an optical Kerr gate or a streak camera could be used to detect early-arriving photons and reject delayed diffusive photons [6]. Spatial filtering of diffusive photons with confocal gating has also been proposed to increase gating efficiency [7]. In these implementations, MW experiencing small deflection angles and also BW were detected due to relatively broad time-gating window (10 ps). Therefore, the spatial resolution of the reconstructed images was far from the optical diffraction limit. In a quasi-diffusive regime where the thickness of the scattering medium is smaller than 10 times the transport mean free path, there occurs temporal overlaps of MW with BW. Those MWs having the same flight time as BW also pass through the time-gating window and contribute to a coherent noise. The MW whose arrival time is within the time-gating window (i.e. the duration of reference pulse) set around the flight time of the BW also interferes with the reference beam and contributes to a coherent noise. Various imaging techniques such as optical coherence tomography (OCT) and confocal microscopy have been developed for detecting early arriving waves with an improved time-gating resolution (100 fs or shorter). The coherence property of light wave was used to generate temporal and spatial coherence gating that discriminates BW from all the other scattered waves, thereby realizing the diffraction-limited imaging [8–10].

In the reflection geometry, BW is manifested as a single-scattered wave (SW) that is scattered only once by a target object, but not at all by the scattering medium [11]. It travels straight through the scattering medium in which the object is embedded on the way to and from the object. Unlike in the transmission geometry, SW has specific time delay set by the depth of target and the average refractive index of the scattering medium, and there exists MW traveling at shallower depth than the target and thus arriving earlier than SW. Similar to the case of transmission geometry, various deep-tissue imaging methods have been developed for the selective detection of SW by using either confocal or coherence gate [9,12]. More often than not, both were combined to maximize the gating efficiency [9,11–15]. Continuous efforts have been made to devise additional gating methods independent of existing modalities such as the space gating method based on detecting ultrasound-modulated BW [16].

So far, most experimental studies have focused on the selective filtering of BW (or SW that does not experience any scattering event at the scattering medium enclosing the object of interest) using various gating methods or the temporal dynamics of diffusive photons (i.e. MW in diffusive regime). Little attention has been paid to investigate the spatio-temporal behavior of the BW and MW at the same time. In fact, the relative temporal behavior of BW and MW remains unclear in the quasi-diffusive regime where the BW is superimposed in the strong MW background in both temporal and spatial domain, albeit its importance in deep tissue imaging. Here, we propose a method to separately quantify BW and MW as a function of flight time and the thickness of the scattering medium. We implemented a low-coherence interferometric microscope that can detect a time-gated transmission matrix of an object embedded within scattering medium [17,18]. Using the measured transmission matrix, we obtained the coherently accumulated BW images by applying the momentum conservation property of the SW reported earlier in the reflection geometry [11] as well as the conventional incoherent image. By comparing them at each flight time, we could separately map out BW and MW even when MW is much stronger than BW. Based on these results, we measured the flight time difference of BW with respect to the most probable flight time of MW and the intensity ratio between BW and MW as a function of the thickness of the scattering medium. The improved understanding of the spatio-temporal behavior of BW and MW will facilitate the deep-tissue and high-resolution imaging of targets under a highly scattering medium.

2. Experimental setup

The experimental setup is based on an off-axis low-coherence interferometer that is designed for recording the complex field of the transmitted wave over a wide field of view (Fig. 1). A mode-locked Ti:Sapphire laser (center wavelength $\lambda = $ 800 nm, bandwidth: 25 nm, coherence length: 30 $\mathrm{\mu }\textrm{m}$) is split into a sample beam and a reference beam by a beamsplitter (BS1). The sample beam was steered by a galvanometer mirror (GM) located at the plane conjugate to the object embedded within a scattering sample and delivered to the sample via an objective lens (OL_ill, Olympus RMS20X, 0.4NA). The sample beam transmitted through a sample was then collected by another objective lens (OL_col, Olympus RMS20X, 0.4NA) and relayed to the camera (sCMOS camera, PCO edge 4.2) located at the plane conjugate to the object. The magnification from the object to the camera plane was 20 and the recorded field of view was $W \times W = 65 \times 65\; \mathrm{\mu}{\textrm{m}^2}$. The reference beam was introduced to a diffraction grating (80 grooves/mm), and the first-order diffracted beam was combined with the sample beam at the camera to generate an off-axis interference image. In the reference arm, a reference beam pathlength scanning mirror (RM) was installed on a translation stage to control the optical path length of the reference beam relative to the path length of the sample beam. Interferograms were recorded at N different incidence angels at the speed of 100 frames per second (fps).

 

Fig. 1. Schematic of the experimental setup. Light source: custom-built Ti:Sapphire laser (center wavelength: 800 nm, bandwidth: 25 nm). GM: galvanometer scanning mirror, BS1-3: beam splitters, RM: reference scanning mirror, G: transmission grating, OL_ill: illumination objective lens, OL_col: collection objective lens, Camera (sCMOS camera). The sketches of temporal intensity profiles are shown for the laser output ${I_{\textrm{in}}}$, sample beam arriving at the camera ${I_\textrm{S}}$, and the reference beam ${I_\textrm{R}}$.

Download Full Size | PDF

The complex field of the sample wave at the time delay $\tau $ for each incidence angle (i.e. each input transverse wave vector ${{\mathbf k}_{\textrm{in}}}$), $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$, was then extracted by performing the Hilbert transform on the measured interferogram [19]. Here ${\mathbf r} = ({x,y} )$ is the spatial coordinates of the sample plane, and $\tau $ is the flight time of the sample wave relative to the reference wave which is controlled by the position of RM. ${{\mathbf k}_{\textrm{in}}} = ({2\pi /\lambda } )\textrm{sin}{{\mathbf \theta }_{\textrm{in}}}$ is the input transverse wavevector for the incidence angle ${{\mathbf \theta }_{\textrm{in}}} = ({{\theta_x},{\theta_y}} )$ at the sample plane. We constructed a time-gated transmission matrix $t({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ for each time delay $\tau $ by assigning $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ to the column and row indices corresponding to ${{\mathbf k}_{\textrm{in}}}$ and , respectively. For the complete coverage of the free propagating modes on the input sides, we scanned ${{\mathbf \theta }_{\textrm{in}}}$ in such a way that the corresponding wave vectors ${{\mathbf k}_{\textrm{in}}}$ covers the aperture of $\alpha $ on a Fourier plane (i.e. $|{{\mathbf k}_{\textrm{in}}}|\le \alpha ({2\pi /\lambda } )$)) with a step size less than $2\pi /W$. Typically, we chose $\alpha = 0.35$, which sets the number of incident angles $N = 2500$.

For each measured time-gated transmission matrix, we applied the collective accumulation of single scattering (CASS) algorithm. CASS was previously developed and employed in the reflection geometry to enrich the ballistic wave contribution to the transmitted waves [11]. The key concept is to first take the Fourier transform of $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ with respect to ${\mathbf r}$ for obtaining the spatial frequency spectrum $\tilde{E}({{\mathbf k},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ and then to perform the coherent summation for various ${{\mathbf k}_{\textrm{in}}}$ in the momentum difference plane. Mathematically, the measured transmitted wave in the frequency domain is given as the summation of the ballistic and multiple-scattered waves:

Here, the ballistic wave can be described as ${\tilde{E}_\textrm{B}}({{\mathbf k},\tau ;{{\mathbf k}_{\textrm{in}}}} )= \sqrt {\mu (L )} \tilde{O}({{\mathbf k} - {{\mathbf k}_{\textrm{in}}}} )$, where $\tilde{O}({\mathbf k} )$ is the spectrum of an object function $O({\mathbf r} )$, the spatial amplitude transmittance of an object at ${\mathbf r}$. $\mu (L )= {\textrm{e}^{ - L/{l_\textrm{s}}}}$ is the intensity attenuation of the ballistic wave traveling through scattering medium with thickness L with the scattering mean free path ${l_s}$. Therefore, the ballistic waves that are singly interacted with the object are invariant in the momentum difference plane. In contrast, the multiple-scattered waves can be considered as uncorrelated speckle fields because the multiple scattering paths for different ${{\mathbf k}_{\textrm{in}}}$ are uncorrelated.

Then, we can rewrite Eq. (1) as a function of the momentum difference $\mathrm{\Delta }{\mathbf k} = {\mathbf k} - {{\mathbf k}_{\textrm{in}}}$ as following

Note that the ballistic contribution is increased by N times in its amplitude because ${\tilde{E}_\textrm{B}}({\mathrm{\Delta }{\mathbf k} + {{\mathbf k}_{\textrm{in}}},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ is independent of ${{\mathbf k}_{\textrm{in}}}$. That is, $\tilde{E}({\mathrm{\Delta }{\mathbf k},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ is coherently added for different ${{\mathbf k}_{\textrm{in}}}$, while the multiple-scattered waves ${\tilde{E}_\textrm{M}}({\mathrm{\Delta }{\mathbf k} + {{\mathbf k}_{\textrm{in}}},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ are incoherently added. Taking the inverse Fourier transform of ${\tilde{E}_{\textrm{CASS}}}\left( {\mathrm{\Delta }{\mathbf k},\tau ;{{\mathbf k}_{\textrm{in}}}} \right)$ with respect to $\mathrm{\Delta }{\mathbf k}$, we obtained ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$. The intensity of the CASS image, ${I_{\textrm{CASS}}}\left( {{\mathbf r},\tau } \right) = {\left| {{E_{\textrm{CASS}}}\left( {{\mathbf r},\tau } \right)} \right|^2}$, is approximately given as $\langle {I_{\textrm{CASS}}}\left( {{\mathbf r},\tau } \right) \rangle \simeq {N^2}{I_\textrm{B}}\left( {{\mathbf r},\tau } \right) + N{I_\textrm{M}}\left( {{\mathbf r},\tau } \right)$, where ${I_\textrm{B}}({{\mathbf r},\tau } )= \langle {|{{E_\textrm{B}}({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )} |^2}\rangle_{{{\mathbf k}_{\textrm{in}}}} = \alpha (L ){|{O({\mathbf r} )} |^2}$ is average intensity of each BW image and ${I_M}({{\mathbf r},\tau } )= \langle{|{{E_M}({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )} |^2}\rangle_{{{\mathbf k}_{\textrm{in}}}}$ is that of the MW [20]. Here $\langle \rangle_{{{\mathbf k}_{\textrm{in}}}}$ indicates the average with respect to ${{\mathbf k}_{\textrm{in}}}$. We also obtained the incoherent image

Finally, from $\langle {I_{\textrm{CASS}}}({{\mathbf r},\tau } ) \rangle$ and $\langle {I_{\textrm{incoh}}}({{\mathbf r},\tau } ) \rangle$, we could extract the intensity of the BM and MW waves, ${I_\textrm{B}}$ and ${I_\textrm{M}}$, as follows.

In Eqs. (5) and (6), we assume that the contribution of MW in CASS image is negligible (i.e. ${N^2}{I_\textrm{B}} \gg N{I_\textrm{M}}$), which is valid when the object is clearly visible in CASS image. Otherwise, it is difficult to discuss ${I_\textrm{B}}$ from CASS image since the fluctuation of MW in CASS image is proportional to $N{I_\textrm{M}}$ because of Gaussian random property of speckle fields. Finally, we can define the ballistic to multiple scattering intensity ratio as $\textrm{BMR} = {I_\textrm{B}}/{I_\textrm{M}}$, and the validity condition of our analysis becomes $\textrm{BMR} \gg 1/N$ (i.e. ${I_\textrm{B}}$>>$({1/N} )\times \; {I_\textrm{M}}$).

3. Temporal distribution of ballistic and multiple-scattered waves

The optical path length of MW is expected to be broadened and delayed with respect to that of BW (shown in Fig. 2(a)). Using the framework presented in section 2, we quantitatively investigated this behavior with the time-gated transmission matrix measured over a large range of delay time $\tau $. For this purpose, USAF target embedded in a thick scattering medium was prepared by dispersing 1.0-µm-diameter polystyrene beads in a slab of polydimethylsiloxane (PDMS) layer with 2.5 g/100 ml concentration. The scattering mean free path of this scattering medium was measured to be ${l_\textrm{s}}$ = 36.9 µm. Also, the transport mean free path of this scattering medium was estimated to be ${l_\textrm{t}}$ = 461 µm. The thickness of this scattering medium is about 2.2 time thicker than the ${l_\textrm{t}}$, which is typically defined as quasi-diffusive regime where the BW is temporally overlapped with the MW. We controlled the magnitude of MW by changing the thickness of the scattering medium.

 

Fig. 2. Time-gated BW and MW images extracted from time-gated transmission matrices. (a) Illustration of the behavior of BW and MW transmitted through a scattering medium. MW is elongated in its flight time distribution. (b) Average BW intensity ${I_\textrm{B}}$ (solid red curve) and MW intensity ${I_\textrm{M}}$ (solid blue curve) detected at camera plane. Both are normalized by the peak of multiple-scattering intensity. For better visualization, ${I_\textrm{B}}$ was multiplied by 177 (red dashed curve). Thickness of scattering medium is 1000 $\mathrm{\mu}\textrm{m}$. (c) Intensity maps of BW and MW as a function of the relative flight time $\mathrm{\Delta }\tau = \tau - {\tau _\textrm{B}}$. ${I_\textrm{B}}$ and ${I_\textrm{M}}$ for each $\mathrm{\Delta }\tau $ were normalized by the maximum intensity in ${I_\textrm{M}}$. Scale bar: 10 µm.

Download Full Size | PDF

In Fig. 2(b), we present the behavior of the average BW intensity ${I_\textrm{B}}$ (solid red curve) and MW intensity ${I_\textrm{M}}$ (solid blue curve) detected at camera plane as a function of flight time when the USAF target is embedded in 1-mm-thick scattering medium corresponding to $L = 27{l_\textrm{s}}$. We found that BW signal has a maximum value when the reference beam arrival time was delayed by ${\tau _\textrm{B}}$ = 1.35 ps with respect to the arrival time without scattering medium. This agrees well with the increased optical path length by 1-mm-thick PDMS layer estimated by ${\tau _B} = \mathrm{\Delta }nL/c$ with $\mathrm{\Delta }n = 0.41$, the average refractive index difference between the scattering medium and air, and $L = 1\; \textrm{mm}$. The peak of MW intensity was delayed by 80 fs with respect to ${\tau _\textrm{B}}$. Also, the peak intensity of MW was 177 times stronger than that of BW, resulting in $\textrm{BMR} = 0.006$, which means the BW is obscured by MW in a single interferogram. However, this BMR value was sufficiently large to reconstruct CASS image and estimate ${I_\textrm{B}}$ because the coherent accumulative nature of CASS algorithm. In OCT, such MW generates the well-known speckle noise, which is the main image-degrading factor. The main scope of our paper is to separately identify the contribution of BW and the coherent MW arriving at the same time-gating window by applying additional spatial coherence gating to the measured time-gated transmission matrix based on the CASS algorithm.

Figure 2(c) shows the map of ${I_\textrm{M}}$ and ${I_\textrm{B}}$ for various relative time delay $\mathrm{\Delta }\tau = \tau - {\tau _\textrm{B}}$. We found that the maps of ${I_\textrm{B}}$ clearly reveals the shape of USAF target at $\mathrm{\Delta }\tau = 0$ while the maps of ${I_\textrm{M}}$ show the uniform MW background. With the increase of $\mathrm{\Delta }\tau $, ${I_\textrm{B}}$ was decreased because the ballistic wave has the same temporal character with the incident wave of a finite pulse width, and ${I_\textrm{M}}$ quickly dominates ${I_\textrm{B}}$. Therefore, the contrast of the ballistic image was quickly reduced, and resolving power was lost. To visualize the effect of temporal coherence gating, we also present the maps of ${I_\textrm{B}}$ and ${I_\textrm{M}}$ under continuous-wave (CW) illumination (the last column in Fig. 2(c)) by turning off the mode-locking of the laser system. For the CW illumination, we could not observe any features of the USAF target from the map of ${I_\textrm{B}}$ as MW here contributes in a time-integrated fashion so that $\textrm{BMR}$ is further reduced to the level where the ballistic wave could be clearly identified using CASS algorithm.

4. Finding ballistic waves using wave correlation

In Fig. 2, we extracted BW and MW contributions based on the image intensities of ${I_{\textrm{CASS}}}$ and ${I_{\textrm{incoh}}}$. In this section, we present a method to extract each contribution from the wave correlation with respect to the CASS image at the peak arrival time of BW, ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$. Since ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$ is close to the ideal BW containing the target image, the normalized cross correlation reports the fraction of BW constituting the field that is correlated with ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$. We first applied this analysis to a relatively thin scattering medium of $L = 200\; \mathrm{\mu}\textrm{m}$, which corresponds to $6{l_\textrm{s}}$, and ${\tau _\textrm{B}} = 0.3\; \textrm{ps}$. The red curve in Fig. 3(a) shows the normalized correlation between ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ and ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$, visualizing the contribution of BW in the CASS image as a function of flight time. Since ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ is composed mainly of BW, FWHM of the curve (131 fs) was close to the coherence time of the light source (100 fs). This wave correlation method can be applied to individual complex-field images $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ taken for each illumination angles. The normalized correlation between $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ and ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$ are shown for a few representative illumination angles in Fig. 3(a). Since the correlation values were too small to be compared with that of ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$, they were separately plotted in Fig. 3(b). Considering that the correlation value was around 0.01, BW accounted for only 0.01% of the intensity in each $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$. Even if BW contribution was weak, the wave correlation method could find its flight time distribution in each image. It was observed that smaller illumination angles contain larger BW contribution, and the flight time distribution for large illumination angle was delayed due to the increase of effective optical path length through scattering medium.

 

Fig. 3. Finding ballistic waves using wave correlation with CASS image at $\tau = {\tau _\textrm{B}}$. (a) Normalized cross-correlation between ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ and ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$ (red line), and normalized cross-correlation between $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ and ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$ for each ${{\mathbf \theta }_{\textrm{in}}}$ for the scattering medium with thickness $L = 200\; \mathrm{\mu}\textrm{m}$ and . (b) Same as (a), but only for to better visualize the individual angular image’s correlation. (c) and (d) Same as (a) and (b), respectively, but when $L = 1000\; \mathrm{\mu}\textrm{m}$ and ${\tau _\textrm{B}} = 1.35\; \textrm{ps}$.

Download Full Size | PDF

We repeated the same analysis for a scattering medium with the increased thickness ($L = 1000\; \mathrm{\mu}\textrm{m}$, $27{l_\textrm{s}}$). As shown in Fig. 3(c) and (d), flight time distribution of BW could hardly be found in each $E({{\mathbf r},\tau ;{{\mathbf k}_{\textrm{in}}}} )$ due to the significant increase of MW. However, BW contribution could still be identified in the correlation between ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ and ${E_{\textrm{CASS}}}({{\mathbf r},{\tau_\textrm{B}}} )$ (red curve in Fig. 3(c)). This proves that the preferential increase of BW contribution in ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ led to the better identification of BW. It is noteworthy that the correlation curve was slightly different from that for thinner scattering medium, and its FWHM (120 fs) was measured to be slightly smaller than that of the thinner scattering medium. This is because the relative contribution of MW was not negligible even in ${E_{\textrm{CASS}}}({{\mathbf r},\tau } )$ as $\tau \; $was set farther away from ${\tau _\textrm{B}}$. In other words, when $|{\tau - {\tau_\textrm{B}}} |> \sim 50\; \textrm{fs}$ in the case of 27${l_\textrm{s}}$-thick sample, the validity condition, ${I_\textrm{B}}/{I_\textrm{M}} \gg 1/N$, of the CASS algorithm is gradually compromised. Therefore, CASS image reconstruction fidelity is decreased, leading to the underestimation of the BW intensity.

5. Relative temporal delay between ballistic and multiple-scattered waves

We investigated the temporal behavior of BW and MW with the increase of scattering medium thickness L (Fig. 4). Figure 4(a) shows the difference between the most probable flight time of MW, the flight time at which MW intensity is maximum, and the BW flight time ${\tau _\textrm{B}}$. With the increase of $L$, this time difference grew drastically, suggesting that the temporal coherence gating is getting more efficient for thicker scattering medium. BMR, the intensity ratio between BW and MW, at ${\tau _\textrm{B}}$ is plotted in Fig. 4(b) for the cases with and without coherence gating. BMR was decreased with the increase of $L$ regardless of temporal gating. In fact, the fidelity of the BMR estimation was low in the case without coherence gating beyond the thickness of $L = 25{l_\textrm{s}}$. This is because the CASS imaging failed to work without coherence gaging Fig. 1(c); right), and the estimation of single scattering intensity is not valid in this case. On the contrary, BMR estimation was reliable all the way to $L = 27{l_\textrm{s}}$ in the case with the time gating. Nevertheless, we could directly estimate the effect of temporal gating from the temporal profiles of BW and MW shown in Fig. 2(b). We integrated MW intensity profile (blue curve in Fig. 2(b)) over all the flight time to obtain the MW detected without time gating. We also integrated MW intensity profile overlapped with the temporal gating window provided by the BW intensity profile (red dashed curve in Fig. 2(b)) to obtain the time-gated MW. Their ratio was the direct measure of the effectiveness of the temporal gating, and it was measured to be 2.2.

 

Fig. 4. Relative temporal delay between BW and MW, and BMR depending on the scattering medium thickness. (a) Time difference between ${\tau _\textrm{B}}$ and the most probable time delay of MW, as a function of $L/{l_\textrm{s}}$. (b) BMR, the ratio between BW intensity and MW intensity, at ${\tau _\textrm{B}}$ depending on $L/{l_\textrm{s}}$. Dashed line indicates the threshold sensitivity level above which the BMR estimation is reliable in the case with no temporal gating.

Download Full Size | PDF

6. Discussion

In this study, we designed an optical system to measure a time-gated transmission matrix in the transmission geometry. The complex-field maps of the transmitted waves were measured for various illumination angles, and they were used to construct the time-gated transmission matrix. We applied the collective accumulation of single scattering algorithm to the transmission matrix measured for each flight time to obtain an object image with the weighting of ballistic wave over multiple-scattered wave by a factor of the number of incidence angles. In addition, we incoherently added all the incidence angle-dependent images to obtain an image that are equally weighted for both ballistic and multiple-scattered waves. By comparing the two images, we could separately obtain the ballistic wave intensity and multiple scattering intensity at each flight time. In doing so, we could map the temporal distributions of ballistic and multiple-scattered waves. Furthermore, we presented the wave correlation method to find the ballistic wave contribution to the transmitted wave taken for each incidence angle. We investigated the temporal characteristics of ballistic and multiple-scattered waves for various thickness of scattering medium ranging up to 27 times the scattering mean free path. The detailed analysis was carried out to find the flight-time difference between ballistic and multiple-scattered waves and single- to multiple-scattering intensity ratio depending on the thickness of the scattering medium. This provided the quantitative understanding of the wave propagation through the scattering medium.

Our presented method is, to our best knowledge, the first method that enables direct observation of temporal characteristic of ballistic waves and multiply scattered waves in a quasi-diffusive regime. Traditional method to detect the ballistic wave was to illuminate a sample with a collimated beam and collect the beam through two apertures located at a large distance [12]. This approach often fails when the thickness of the medium ($L/{l_\textrm{s}}$) becomes larger than 10. Although small in magnitude, the multiply scattered wave that passes through the apertures becomes larger than the exponentially decaying ballistic wave. To circumvent this issue, time gating approach has been devised, but this approach is only successful in a fully diffusive regime where the difference in the arrival time of the ballistic wave and multiply scattered wave is sufficiently larger than the temporal resolution of the time gating. Therefore, in quasi-diffusive regime where the ballistic wave temporally and spatially overlaps with the multiply scattered wave, there has been difficulties in discriminating ballistic and multiply scattered wave. We overcome this challenge by coherently accumulating, thereby selectively amplifying, the ballistic wave buried within an intense scattered wave. The amplification factor is given by the number of orthogonal free modes incident on the sample, which is closely related with the key parameters of our measurement system – field of view and numerical aperture. The number of free modes is linearly increased with the field of view or the aperture area on a Fourier space (i.e. the square of NA).

In our previous study, we implemented the coherence gating and coherent accumulation in a reflection geometry [11]. It is worth noting that there are slight differences in the role and implementation of two operations - coherence gating and coherent accumulation - in the transmission and reflection geometries due to a geometrical reason even though the coherence gating plays the role in rejecting the multiple-scattered waves and the coherent accumulation then provides the depth selectivity along with further suppression of multiple-scattered wave based on spatial coherence in both the transmission and reflection geometry. In reflection geometry, the reference arm length is set to the midst of temporal profile of arriving waves and is dependent on the object position. The coherence gating in reflection geometry not only rejects the multiple-scattered waves as in the transmission geometry, but also provides depth selectivity on the single-scattered waves. On the contrary, the role of coherence gating in transmission geometry is to solely reject the multiple-scattered waves. For this purpose, the reference arm length is set to the earliest flight time of the arriving waves regardless of the object position. Then coherent accumulation exclusively provides depth selectivity.

The newly presented wave correlation method for finding the ballistic waves has similarity with the traditional two-aperture method described above. As if the second aperture filters the collimated input defined by the first aperture, wave correlation with respect to the CASS image selectively detects ballistic waves preserving the image information. However, the presented method is more general in the sense that it is applicable to analyzing the ballistic waves associated with the image formation. Furthermore, the use of CASS imaging amplifying the ballistic waves raises the sensitivity of the finding the ballistic wave, thereby extending the working range well beyond 10 $L/{l_\textrm{s}}$. The effect of temporal gating is rather mild in the present study. Setting $\langle \theta\rangle $ as the average deflection angle of the scattered photons propagating inside a medium to the incident angle, it is often small in quasi-diffusive regime. Then, the path length difference relative to that of ballistic photons, given as $L\left( {\frac{1}{{\cos \langle\theta\rangle }} - 1} \right)$, can be approximated to $L{\langle\theta\rangle ^{2}}/2$. Even with $\langle\theta\rangle $ up to 0.3 rad, the path length is expected to be ∼ 45 $\mathrm{\mu}\textrm{m}$ for $L = 1\; \textrm{mm}$. In our experiment, the effect of temporal gating was not very significant because this difference is comparable to the coherence length (30 $\mathrm{\mu}\textrm{m}$) of our system. However, one may expect that the effect of time gating defined by the intensity ratio of BW and MW will be dramatically increased for either case where L or anisotropy factor g is increased, or the temporal gating window becomes narrower. Conversely, we may conclude that the multiply scattered wave that are captured near the flight time of the ballistic wave are mostly comprised of the waves with small deflection angles, often called snake photons. Varying time-gating window may provide us with a route to noninvasively deciphering the deflection angles of multiply scattered waves within the scattering medium. Our method, as a first study to differentiate the ballistic and multiply scattered waves in scattering medium in quasi-diffusive regime, will pave a way toward developing new imaging modalities that deterministically utilize the snake photons.

In a longer term, our study will help understanding the behavior of the snake waves that have a similar flight time with the ballistic wave as they propagate the scattering medium in a forward scattering trajectories. We expect the better understanding will result into the new methodologies for deeper imaging that overcomes the limit of the conventional methods that only rely on the exponentially decaying ballistic wave. Finally, such new imaging methods will enable thick-tissue pathology. Conventional biopsy techniques rely on the complicated and time-consuming preprocessing of extracted tissues including freezing, sectioning, fixation, and staining. In the long run, our coherent imaging method can visualize important pathological features within thick fresh tissues without labeling and other processes, making it a vital tool for pathology and the associated applications in medical diagnosis.