Coherent fluctuations in time-domain diffuse optics

Optical waves multiply scattered by biological tissues contain useful information about the sample constituents and their microstructures.¹ Earlier, it was discovered that light in the near infrared (NIR) region can probe the tissue down to 2 cm–3 cm in depth. In fact, in that spectral window, absorption is much smaller than scattering; thus, photon propagation can be efficiently described as a diffusion process. Many diffuse optical methods were then developed for probing tissues in a non-invasive way.^2,3

When a narrow-band light source is used, photons emerging from a random medium acquire path-length-dependent phases that cause strong spatial variations of intensity, or a so-called speckle pattern.^4,5 Even if light—or waves in general—are multiply scattered, speckle effects are observed and exploited in many fields, allowing imaging^6–8 or monitoring⁹ the quantities not directly observable. If scattering elements move, the speckle intensity fluctuates over time due to slight changes in photon paths. Diffuse correlation spectroscopy (DCS) is an optical technique that exploits this effect, which, using a continuous-wave laser source, measures the temporal intensity decorrelation related to particle movements.¹⁰ In biomedical optics, it finds an important application to measure the motion of particles such as red blood cells, and thus, blood flow (BF).¹¹ 

With the use of long-coherence-length pulsed lasers, recently, it was possible to record intensity decorrelation at different photon times-of-flight. Thanks to the physical relationship between the time-of-flight and mean-penetration-depth,¹² this technique, called time-domain (TD) DCS, demonstrated the discrimination of blood flow at different depths, both in small animals¹³ and in humans.¹⁴ A different approach for resolving scatterer dynamics at different times-of-flight is the use of an interferometric setup together with a wavelength-swept laser, being the basis of a recent technique called interferometric NIRS (iNIRS).¹⁵ The assumption behind TD-DCS is that the single-photon measurement freezes a microscopic state faster than the sample dynamics. For this reason, it is possible to define two independent temporal scales: the photon time-of-flight so-called micro-time, corresponding to the light propagation time scale (∼ps to ns), and the elapsed time or macro-time, corresponding to the particle dynamics time scale (∼μs to ms). This permits us to interpret time-of-flight localized intensity-changes as localized particle movement similarly to acoustic wave seismology,⁹ giving spatial information on the dynamics, e.g., the Brownian diffusion coefficient of the scatterers. Thus, in biological tissues, it may give additional information regarding the micro-vascular blood flow.¹¹ 

In this work, we report the observation of temporal fluctuations encountered within a single distribution of times-of-flight (DTOF) curve, deviating from the behavior predicted by the diffusion theory. This effect becomes evident for a high temporal coherence laser source, and for samples composed of a stiff structure. We will show that the DTOF fluctuations are well described by speckle statistics¹⁶ and their evolution with respect to the elapsed time (macro-time) follows a fluctuating behavior. For this reason, this phenomenon may be used not only to extract depth-resolved information regarding the sample dynamical properties but also to evaluate the coherence properties of the light source. Understanding this phenomenon is a crucial task for the correct estimation of the sample optical properties, i.e., scattering and absorption coefficients, from the measured DTOF.³ The combined assessment of static (absorption and scattering) and dynamic (blood flow) properties provides, for instance, a more complete description of hemodynamics in the brain or in cancerous lesions.¹ 

Figure 1(a) shows a scheme of the experimental TD-DCS setup used for the measurement, described in more detail in Ref. 14. It is based on a custom-made active mode-locked Ti:Sapphire (Ti:Sa) laser with 100 MHz repetition rate, tuned to the wavelength λ = 785 nm, and pulse full-width at half-maximum (FWHM) of 200 ps. Light with 50 mW average power is delivered to the sample with a graded index fiber, recollected in reflectance geometry at the source–detector separation ρ = 1.5 cm with a single-mode fiber, and sent to a 50 µm diameter single-photon avalanche diode (SPAD). Photon arrival times, so-called time-stamps, are then acquired by a time-correlated single photon counting board (TCPSC) with 25 ps temporal resolution. Then, from the time-stamps, we compute the distribution of times-of-flight (DTOF) with an integration time of 100 ms, representing the binning time of our measurements in the macro-time direction. The system is characterized by measuring the instrument response function [IRF, black curve in Fig. 1(b)], acquired by facing the injection and detection fibers separated by a thin (<0.5 mm) Teflon sheet.

With this experimental setup, we measured the solid resin sample described in Ref. 17 (4.5 cm height and 10.5 cm diameter) designed to have nominal reduced scattering coefficient μ′_s = 18.6 cm⁻¹ and absorption coefficient μ_a = 0.009 cm⁻¹. Figure 1(b) reports the measured DTOF averaged over 100 repetitions, i.e., 10 s (red curve), which exhibits strong fluctuations with respect to the time-of-flight. Under the diffusion approximation, the theoretical reflectance can be expressed as

where k is a proportionality constant, v is the speed of light in the medium, and S is a boundary condition factor, which can be calculated using the method of images.¹⁸ 

Longer averaging times (5 min) strongly reduce the fluctuations and permits us to fit the DTOF with the theoretical time-resolved reflectance convoluted with the IRF (blue curve), with respect to optical properties and an arbitrary time shift.^17,19 The long integration time is needed to effectively average the detected intensity over several microscopic configurations for comparing the measurement with the diffusion theory that neglects any interference effect in the detected fields, as demonstrated for acoustic waves.²⁰ Performing the fit in the region between 80% on rising and 5% on falling edge of the DTOF, the retrieved optical properties are $μs′$ = 21.2 cm⁻¹ and μ_a = 0.010 cm⁻¹ similar to the expected values.

To elucidate the physical origin of the speckle deviation with respect to the diffusion theory, we studied its dependence on the sample rigidity considering three different phantoms: (a) a 5:95 liquid mixture of lipid (Intralipid 20%) in distilled water, nominal properties μ′_s = 10 cm⁻¹ and μ_a = 0.02 cm⁻¹, (b) a flexible silicone phantom, having μ′_s = 10 cm⁻¹ and μ_a = 0.05, cm⁻¹, and (c) the previous solid resin phantom. For each sample, 3000 reflectance curves were acquired with a sampling time of 100 ms. In Figs. 2(a)–2(c), we report for the three mentioned samples, the relative variation (i.e., contrast) of the measured reflectance (R) with respect to its temporal average (⁠$⋯τ$⁠), defined as

where t is the time-of-flight (i.e., micro-time, vertical axis in Fig. 2) and τ is the elapsed time (i.e., macro-time, horizontal axis in Fig. 2). We studied the behavior of this fluctuation also as a function of the detection fiber’s core diameter. In particular, Figs. 2(d)–2(f) report the measured contrast C(t, τ) for another solid resin phantom with optical properties μ′_s = 10 cm⁻¹ and μ_a = 0.1 cm⁻¹ with decreasing values of the core diameter: 25 μm, 10 μm, and 5 µm.

As shown in Fig. 2 (top row), the fluctuations contrast strongly increases when samples with slower dynamics are considered, in particular moving from fluid [Fig. 2(a)] to solid samples with increasing stiffness [Figs. 2(b) and 2(c)]. In addition, Fig. 2 (bottom row) illustrates the dependence of the measured contrast with respect to the detection fiber diameter, thus on the number of speckles detected simultaneously. Increasing the collection area implies averaging through a higher number of spatial speckle modes [Fig. 2(d)], reducing the visibility of the fluctuation.^21,22 Their contrast strongly increases when the detection fiber diameter approaches the average speckle size [Figs. 2(e) and 2(f)].

In order to verify quantitatively the hypothesis that the observed fluctuations are an effect of temporal speckle, we apply the model developed by Goodman, who analyzed the photo-counts probability distribution in the presence of speckle and Possion noise.^16,23 The time-resolved variance $σI2(t)$ of the measured DTOF may then be written as

In Eq. (3), the time-resolved variance $σI2(t)$ is given by the sum of two contributions: Possion noise (first term), proportional to the average intensity, and speckle fluctuations (second term), proportional to the square of the average intensity. The coefficient β(t), often called coherence parameter, is inversely proportional to the detected effective number of modes, which depend on the detection fiber’s design (e.g., core diameter, numerical aperture) and the collected wavelength.²¹ It is related to the amplitude of the intensity fluctuations and it can be interpreted as the amount of coherence, both temporal and spatial, in the detected fields.²² For non-polarized detection, its maximum theoretical value is 0.5 since the speckle patterns of the two independent polarizations add up incoherently, halving the resulting intensity fluctuations.⁴ 

Regarding Eq. (3), we stress that it is rigorously valid only in the quasi-static regime, i.e., under the assumption that the scatterer dynamics is slow compared to the measurement integration time since it does not take into account the decrease of speckle contrast due to the sample de-correlation. For this reason, in the following, we will apply it for testing the fluctuation statistics in the case of solid phantoms.

We evaluate the fluctuation statistics of the time-resolved intensity by measuring a solid resin phantom with nominal optical properties μ′_s = 4.7 cm⁻¹ and μ_a = 0.059 cm⁻¹ tuning the pulse width, quantified by the IRF FWHM, to three different values: 450 ps, 225 ps, and 175 ps. The pulse width was modified, acting on the radio-frequency power of the acousto-optic modulator, in order to study the dependence of the DTOF fluctuations on the source temporal coherence. For each pulse width, we computed the time-resolved intensity auto-correlation functions g₂(t, τ) from the measured time-stamps, using 100 ps-wide temporal gates. The intensity auto-correlation function may be expressed as follows:^13,24

In the first step of Eq. (4), we used the Siegert relation, which relates the intensity auto-correlation function to the electric field auto-correlation g₁(t, τ).²⁵ In the second equality, we used the expression of the time-resolved auto-correlation g₁(t, τ), where k₀ is the light wave-number in the medium, D_B is the Brownian diffusion coefficient of the scatterers, and α is the fraction of moving to total scatterers. Thus, the amplitude of g₂(t, τ) is determined by the β(t) parameter, while its decay-rate by the sample dynamics (i.e., D_B).

For extracting robustly the coherence parameter β(t) from the measured auto-correlations g₂(t, τ), we computed, for each time point t, the amplitude of their best fit with a decreasing exponential function, see Eq. (4), considering τ values between 10⁻⁶ and 10⁻⁵ s. For these values of τ, the auto-correlation function is close to its maximum value, and the effect of the after-pulsing (present at shorter τ) is limited.¹⁴ In Figs. 3(a)–3(c), we compare, for the solid resin phantom, the coherence parameter measured from g₂(t, τ) (blue circles with shaded standard deviation) with the one retrieved using Eq. (3) with the measured variance $σI2(t)$ and average intensity $I(t,τ)τ$ (dashed lines), for all the three considered IRF FWHM. After the computation, for both methods, we performed a 200 ps moving average on the micro-time direction.

For all the three conditions in Fig. 3, the measured β(t) decreases with an increase in the times-of-flight t due to the combined effect of IRF tail and finite coherence time.^14,26 The initial value β(t = 0) is equal to the theoretical limit (0.5) for the broader FWHM, Fig. 3(a), while it decreases to ∼0.2 for the narrower FWHM, Figs. 3(b) and 3(c). Within the experimental variability, the value of β(t), retrieved from Eq. (3) with the measured intensity fluctuations, well matches the value measured directly from the amplitude of the auto-correlations, for all the three laser set-points, which indicates that the observed fluctuations can be described by a mixture of pure speckle statistics with Possion noise.

The discrepancy of the estimated from the measured β(t) around time zero may be due to the scarcity of photons at the beginning of the DTOF rising edge, i.e., when it emerges from the noise background, which biases the variance estimation. However, after the first ∼100 ps there is no statistical difference between the two. We emphasize that in Fig. 3, we compared the intensity statistics to the theory using a solid phantom, since Eq. (3) is valid only in the quasi-static regime. For samples exhibiting dynamical properties, like for instance, liquid phantoms, the reduction of the variance of the intensity due to scatterers motion has been studied.²⁷ 

This work reports the first experimental evidence of coherent fluctuations in the field of time-domain diffuse optics. This effect was observed with a >100 ps mode-locked Ti:Sapphire laser, rigid samples, and a single mode detection fiber [Figs. 1 and 2(c)]. The disappearance of the fluctuations for liquid and, to a lesser extent, flexible samples [Figs. 2(a) and 2(b)] may be explained by the self-averaging effect induced by scatterers’ motion. Still, in the case of fluid samples, the effect may be observable with shorter (<ms) integration times, with the need for an increased system throughput. We also observed that the effect is present only when the number of modes of the detection fiber, and thus, the number of coherence areas detected simultaneously, are decreased to very few units [Figs. 2(d)–2(f)]. In this case, the disappearance of the effect can be explained by the presence of spatial averaging. Collectively, these results give evidence of the coherent nature of the observed fluctuations.

We note that similar experimental evidences of these fluctuations were reported in other fields, such as laser ranging,²³ acoustic waves,²⁰ exciton resonant emission,²² and THz spectroscopy,²⁸ and were sometimes named time-resolved speckles. Similar to the case considered in this Letter, the evidence of temporal speckle in diffusive media was observed as a function of the observation angle^29,30 in single-shot measurements (shorter than the characteristic speckle dynamics), while here, as a function of the elapsed time. This reinforces the conclusion that increasing the numerical aperture of the detection effectively acts as a spatial average over different speckle realizations [as in Figs. 2(d)–2(f)].

We have shown the dependence of the phenomena on the sample dynamical properties, and that the measured photo-count distributions can be well described by speckle statistics (Fig. 3). Thus, the statistical properties of the detected light can also be useful in determining the amount of coherence of the diffused pulse, or for marking deviations from ergodic conditions. We note that the reported experiments were limited to samples with homogeneous dynamics. In addition, similar to classical DCS, the evolution of the speckle pattern may be studied to infer information on the sample scatterer particles. However, the additional dimension of the photon time-of-flight, related to its depth penetration probability,¹² may ease the detection of dynamical changes deeply buried in the diffusive media.

With this contribution, we aimed to give physical insights on the propagation of pulsed waves in random media, which may be useful for biomedical applications and material science. We have shown that, in the experimental conditions considered, the photon diffusion equation describes accurately only the temporal average of the detected fields. We have also provided a statistical method to quantify the time-resolved coherence fraction of the re-emitted light, which may be also exploited in many imaging/tomographic techniques.