
1.IntroductionDiffuse correlation spectroscopy (DCS) is an optical technique that measures the scatterer’s motion in diffusive media, traditionally, by injecting a coherent continuous wave (CW) laser beam in a turbid sample and by characterizing the speckle fluctuations by their intensity autocorrelation function. Afterward, using a correlation diffusion model, the motion of red blood cells in the microvasculature is estimated, which has been shown to be proportional to local blood flow (BF).^{1} The technique is an advancement over diffuse wave spectroscopy,^{2} which is widely used in the study of nonbiological media. Over the last two decades, DCS has emerged as a promising technique for many clinical applications since it provides a direct measure of BF in deep ($>1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$) tissues. In standard DCS, which uses a CW light source, photons that travel all the path lengths in the tissue are detected without discrimination. On the contrary, using a pulsed light source, it is possible to measure the photon timeofflight (ToF) and, thus, to select a specific range of path lengths. By exploiting the relationship between photon ToF and mean penetration depth,^{3}^{,}^{4} depthdiscrimination in the BF can be enabled, which is the first advantage of the newly emerging timedomain (TD) DCS with respect to standard (i.e., CW) DCS.^{5}^{,}^{6} A second advantage is the parallel estimation, using the same probe, of the tissue optical properties. This has been shown to increase the accuracy of the BF estimation^{7} and allows, by using more than one wavelength, the determination of oxy and deoxyhemoglobin concentrations alongside BF. The first proofofconcept of TDDCS was given by Yodh et al.^{8} In this work, by using a gating scheme based on an optical delay line and a nonlinear crystal, they showed that it is possible to measure path length resolved autocorrelation functions. However, their scheme was quite complex and needed a high optical power, unsafe for in vivo applications. Recently, a demonstration on phantoms and small animals was given by Sutin et al.,^{5} where a pulsed laser and a timecorrelated singlephoton counting (TCSPC) module were used to measure the BF with narrow (tens of picoseconds) time gates. This allowed the researchers to overcome an important barrier, namely the limited temporal coherence of a pulsed light source. Later, Pagliazzi et al.^{6} overcame the signaltonoise ratio (SNR) barrier for human in vivo measurements, by using a custommade modelocked pulsed light source and broader (few nanoseconds) temporal gates. A hardware gating acquisition scheme, which uses fastgated singlephoton detectors enabling very short source–detector (SD) separation measurements, was also recently demonstrated.^{9} Recently, a theoretical model for TDDCS in multilayer turbid media has been proposed for improving data analysis in nonhomogeneous biological tissues.^{10} Another method to achieve path length resolved DCS, initiated by the works of Tualle et al.,^{11}^{,}^{12} has been recently shown to be available in vivo on small animals.^{13} While in TDDCS, path lengths are resolved based on their ToF, and interferometric measurements achieve the same goal by sweeping source wavelength in time. This promising technique, which can directly measure the electric field autocorrelation function, has not yet been translated to humans. The current challenges of the technique are both on the technological side, regarding the development of effective source and detection schemes, and on the theoretical side, for the improvement of the physical understanding, necessary for data analysis. In this context, our work focuses on the study of finite temporal resolution effects in TDDCS. These effects were very recently studied using a model that is applicable only to narrow time gates.^{14} Here, we argue, as in our earlier work,^{6} that broad gates are crucial to achieve a good SNR and, motivated by this argument, build a more complete mode. This approach would also permit movement toward an accurate depthresolved BF measurement, and can be useful for driving hardware design. As we show, these effects have an important impact on the accuracy of the BF estimates. We also provide and validate a model for BF estimation in the case of a realistic (i.e., nonideal) TDDCS experiment. In the following, we derive the model, based on the instrument response function (IRF), for describing gated autocorrelation functions (Sec. 2); then, we describe the Monte Carlo (MC) simulator, the experimental system, and the phantom model exploited (Sec. 3); then the validation of the model with a numerical simulation (Sec. 4) and a phantom experiment (Sec. 5) are presented. Finally, some discussion and conclusions are reported (Sec. 6). 2.TheoryIn classical CWDCS, for ergodic systems and infinitely coherent sources (i.e., sources with coherence length much greater than the mean path length of photons in the sample), the normalized intensity autocorrelation function ${g}_{2}(\tau )$, where $\tau $ is the correlation delay time, can readily be related to the normalized electric field autocorrelation function ${g}_{1}(\tau )$, by means of the socalled Siegert relation:^{15}^{,}^{16} Here, $\beta $ is a constant, equal to one for an ideal experiment, that is proportional to the inverse of the number of detected modes and decreases for decreasing source coherence length. We note that $\beta $ halved by the use of the typical unpolarized source–detection pairs in DCS experiments. In the diffusive limit, we have^{15}^{,}^{17} where $f(s)=\frac{R(s)}{{\int}_{0}^{\infty}R(s)\mathrm{d}s}$ is the normalized distribution of photon path lengths, often estimated with the theoretical timeresolved diffuse reflectance $R(s)$, and ${g}_{1,s}(\tau ,s)={e}^{ks\tau}$ is the normalized electric field autocorrelation for a single path length $s=vt$, where $v$ is the speed of light in the medium (assumed constant) and $t$ is the photon ToF. Furthermore, the parameter $k({\mathrm{cm}}^{1}\text{\hspace{0.17em}}{\mathrm{s}}^{1})$ is the autocorrelation decay rate per unit path length, which is equal to $2{\mu}_{s}^{\prime}{k}_{0}^{2}\alpha {D}_{\mathrm{B}}$. Here, ${k}_{0}^{2}$ is the square of the wavenumber of light in the medium, $\alpha $ is the fraction of moving scatters, ${\mu}_{s}^{\prime}$ is the reduced scattering coefficient, and ${D}_{\mathrm{B}}$ is the Brownian diffusion coefficient of the scatterers. In most biological media, in fact, photons undergo many scattering events along their path from the source to the detector. Each scattering event with a moving particle gives rise to a fluctuation of the detected speckle intensity, contributing to the decay of the electric field autocorrelation function. Assuming the fields from different paths to be uncorrelated, in Eq. (2), ${g}_{1}(\tau )$ is obtained as a weighted average over the entire path length distribution.In a TDDCS experiment, Eq. (2) is modified to account for the selection of the path lengths in a specific range $[{s}_{0},{s}_{0}+\mathrm{\Delta}s]$, where ${s}_{0}$ is the gate opening path length and $\mathrm{\Delta}s$ is the gate width. The parameter ${s}_{0}$ is chosen to probe the BF at different depths: with earlier or later gate positions, superficial or deeper BF is measured, respectively. Thus, we can write as follows: Eq. (3)$${g}_{1}(\tau )={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}f(s){g}_{1,s}(\tau ,s)\mathrm{d}s.$$In a realistic TDDCS experiment, it is worth noting that the injected light pulse has a finite temporal duration and the detector has a finite temporal response/resolution. When the photon ToFs are recorded (e.g., using a TCPSC module) and are subsequently timegated via software, the system can be assumed to be linearly timeinvariant and, thus, is fully characterized by measuring its IRF. IRF is typically obtained by facing the injection and the collection fibers against each other while taking care to fill all the modes. Notably, the assumption of temporal invariance may not be valid when a fastgated detector is used for gating the photons via hardware, which does so by switching on the detector only in the desired temporal window.^{18}^{,}^{19} We note that Eq. (3) holds only in the case of an ideal system, for which $\mathrm{IRF}=\delta (t)$, where $\delta $ is the Dirac delta function. In the case of an ideal system, it is assumed that it is possible to measure the photon path length $s$ with a perfect accuracy and thus select, by temporal gating of the distribution of timeofflight (DTOF), a specific range $[{s}_{0},{s}_{0}+\mathrm{\Delta}s]$ of path lengths. On the other hand, in a real system, the path length selection is affected by inaccuracies due to the system’s nonideal temporal response. This is because the path length that a photon ideally that undergoes in the medium (i.e., the ideal path length, ${s}^{\prime}$) can differ from the path length that we estimate using the measured ToF (i.e., the measured path length, $s$) in a stochastic manner. In timeresolved spectroscopy,^{19} the spreading of path length is well known and affects the expression of the normalized distribution of path lengths: $\tilde{f}(s)=\tilde{R}(s)/{\int}_{\infty}^{+\infty}\tilde{R}(s)\mathrm{d}s$, where $\tilde{R}(s)={\int}_{\infty}^{+\infty}\mathrm{IRF}({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}$ is the convolution by the IRF of the theoretical timeresolved diffuse reflectance $R(s)$. In our approach, the expression of ${g}_{1,s}(\tau ,s)$, the normalized electric field autocorrelation for a single path length, in the realistic case is calculated as a weighted average of its value in the ideal case. Assuming a fixed measured path length $s$, the number of detected photons within an infinitesimal interval of path lengths between ${s}^{\prime}$ and ${s}^{\prime}+d{s}^{\prime}$ can be expressed as $dN=\mathrm{IRF}({s}^{\prime})R(s{s}^{\prime})d{s}^{\prime}$. We, therefore, have Eq. (4)$${\tilde{g}}_{1,s}(\tau ,s)=\frac{\int {g}_{1,s}(\tau ,s{s}^{\prime})\mathrm{d}N}{\int \mathrm{d}N}=\frac{{\int}_{\infty}^{+\infty}{g}_{1,s}(\tau ,s{s}^{\prime})\mathrm{IRF}({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}}{{\int}_{\infty}^{+\infty}\mathrm{IRF}({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}}.$$Eq. (5)$${g}_{1}(\tau )={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\tilde{f}(s){\tilde{g}}_{1,s}(\tau ,s)\mathrm{d}s.$$Eq. (6)$${g}_{1}(\tau )={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}[{\int}_{0}^{\infty}f({s}^{\prime}){g}_{1,s}(\tau ,{s}^{\prime})\mathrm{IRF}(s{s}^{\prime})\mathrm{d}{s}^{\prime}]\mathrm{d}s.$$Eq. (7)$${g}_{1}(\tau )={\int}_{0}^{\infty}\mathrm{EGF}({s}^{\prime})f({s}^{\prime}){g}_{1,s}(\tau ,{\mathrm{s}}^{\prime})\mathrm{d}{s}^{\prime},$$Eq. (8)$$\mathrm{EGF}({s}^{\prime})={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\mathrm{IRF}(s{\mathrm{s}}^{\prime})\mathrm{d}\mathrm{s}.$$The function $\mathrm{EGF}({s}^{\prime})$ can be interpreted as the probability, for a photon with ideal path length ${s}^{\prime}$, to be mapped into the range $[{s}_{0},{s}_{0}+\mathrm{\Delta}s]$ or equivalently in the corresponding time gate. This function models the fact that, due to finite temporal resolution, also path lengths outside the gate can be assigned to it, contributing to the measured autocorrelation decay rate. To understand the physical meaning of the EGF, Fig. 1 shows two computed EGFs, for an early and a late gate, using an IRF measured with our experimental system.^{6} The measured IRF in Fig. 1 is a composition of the bellshaped excitation pulse and the detector temporal response. The IRF width [full width at halfmaximum (FWHM)] is determined mainly by the input pulse, having a duration of $\sim 300\text{\u2013}400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ps}$ since the time jitter of the detector and the temporal broadening of the gradedindex and singlemode fibers we used are much smaller than the width of the laser pulse. The long diffusion tail is a characteristic of the singlephoton avalanche diode that was used and depends on the detector architecture.^{18} From Fig. 1, it is also possible to see how the ideal gating probability (equal to 1 inside the gate and 0 outside) is modified by the IRF. In fact, the EGF describes, for each time point, how much this ideal gating probability is distorted by the finite temporal resolution of the system. Due to the characteristic tail of the IRF toward positive times, many short path length photons are mapped into the late gate and contribute to its autocorrelation decay rate (see the leading edge of EGF $L$ in Fig. 1). This effect is nonnegligible due to the higher number of early photons compared to late ones.^{3}^{,}^{4} We note that Eqs. (7) and (8) are applicable for any gate width and predict a distortion of the autocorrelation function also in the narrow gate limit used in Ref. 5. In fact, in this limit, by means of the wellknown mean value theorem for integrals, we can write [compared to Eq. (8)] as follows: Eq. (9)$$\mathrm{EGF}({s}^{\prime})={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\mathrm{IRF}(s{s}^{\prime})\mathrm{d}\mathrm{s}=\mathrm{IRF}(\overline{s}{s}^{\prime})\mathrm{\Delta}s,$$Eq. (10)$${g}_{1}(\tau )=\mathrm{\Delta}s{\int}_{0}^{\infty}\mathrm{IRF}(\overline{s}{s}^{\prime})f({s}^{\prime}){g}_{1,s}(\tau ,{s}^{\prime})\mathrm{d}{s}^{\prime},$$3.Materials and Methods3.1.Monte Carlo SimulatorTo study the signals involved in a TDDCS experiment, we have performed MC simulations using software developed at ICFO based on the MCML algorithm.^{20} In particular, we have simulated a semiinfinite homogeneous medium surrounded by air in the reflectance geometry by assigning the reduced scattering coefficient (${\mu}_{s}^{\prime}$), the absorption coefficient (${\mu}_{\mathrm{a}}$), and the refractive index ($n$). We have considered a collimated beam perpendicular to the surface and placed, at a fixed SD separation $\rho $, a 1 mm/side cubic detector. For every detected photon packet $i$, the path length in the medium ${s}_{i}$ and the arrival weight ${W}_{i}$ are recorded in a history file. In postprocessing, the ToF is calculated starting from the path length. To introduce the timing inaccuracy of the system, for every detected packet, we have added to the before mentioned ToF a time shift $\mathrm{\Delta}t$ obtained by performing an MC sampling^{20} of a realistic IRF. First, the IRF cumulative distribution function (CDF) is computed by numerical integration, then the $\mathrm{\Delta}t$, such that the CDF is equal to a given random number $\xi \in [\mathrm{0,1}]$, is selected. Finally, the gated autocorrelation function is computed, in diffusive approximation, in the following way:^{15} Eq. (11)$${g}_{1}(\tau )=\sum _{i=1}^{N}{W}_{i}\text{\hspace{0.17em}}\mathrm{exp}(2{\mu}_{s}^{\prime}{k}_{0}^{2}{D}_{B}{s}_{i}\tau ),$$3.2.Experimental SetupFor experimental validation, we have used the custommade Ti:Sapphire laser source, operated in active modelocking, that was deployed in our previous works.^{6}^{,}^{9} The acoustooptic modulator permits the phaselocking of the longitudinal modes, enabling a pulse repetition frequency of 100 MHz. The wavelength $\lambda $ was tuned to 785 nm. The light was sent to the surface of the sample using a $200\text{}\mu \mathrm{m}$ core stepindex fiber and the diffused light, recollected with a $5\text{}\mu \mathrm{m}$ core singlemode fiber (780 HP, Nufern, East Granby, Connecticut), was delivered to a singlephoton avalanche diode detector (PDM, Micro Photon Devices, Bolzano, Italy). We have implemented optical coupling of the fibers with laser and the detector, respectively, by using a combination of lenses and aperture irises, which were adjusted manually. A synchronization signal was generated using a beamsplitter and photodiode (OCF401, Becker & Hickl, Berlin, Germany) placed at the source’s output. A TCSPC module (PicoHarp 300, PicoQuant, Berlin, Germany) recorded, for every detected photon, two quantities: the photon time stamp, socalled arrival time, and, with a resolution of 8 ps, the delay with respect to the synchronization signal, socalled pulse time. To stay within maximum sync rate of the TCSPC module, which is 84 MHz, a SYNC divider circuit (SDC, developed in our laboratory) was used to reduce the sync rate by two times to 50 MHz. The SDC operates on the photodiode electronic signal by removing the signal generated by exactly one pulse from every two pulses. The resulting TCSPC histogram is composed of two identical replicas, separated by a laser pulse period (10 ns), which are summed together in postprocessing in order not to lose useful photons. More details on the setup can also be found in Ref. 6. The IRF was acquired by facing the input and output fibers, separated by a thin polytetrafluoroethylene sheet. Its FWHM was $\sim 400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ps}$, see Fig. 1. 3.3.Phantom Studies and Data AnalysisWe have mixed a tissuemimicking homogenous liquid phantom by adding 6.5% in mass of Lipofundin (B. Braun Melsungen AG, Germany) to water.^{21} The phantom optical properties were estimated with our system to be ${\mu}_{\mathrm{a}}=0.023\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{1}$ and ${\mu}_{s}^{\prime}=10\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{1}$ by fitting the measured DTOF curve with the theoretical timeresolved diffuse reflectance for a semiinfinite homogeneous medium, convolved with the IRF.^{4} The optical coefficients are then used to compute the normalized path length distribution $f(s)$ by normalizing the theoretical reflectance $R(s)$ to the unitary area. We note that $R(s)$ can also be estimated directly by deconvolving the IRF from the measured DTOF.^{22} This method could be an alternative to estimate the path length distribution, but a high number of counts than the used approach is necessary for a correct inversion to avoid noise in the results. After the experiment, the gated autocorrelation function is computed using the software correlator, as was used in our previous work.^{6} Its working principle can be divided conceptually into two parts. First, the photons with pulse time belonging to the desired temporal window are selected. Then, the arrival times of the selected photons are correlated with the fast photon autocorrelation algorithm described in Ref. 23. 4.SimulationsTo validate the proposed model, we have performed an MC simulation, as described in Sec. 3.1. We have simulated a system with the same optical properties of the phantom described in Sec. 3.3. The simulated values of the other parameters were $\rho =1.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$, $n=1.33$, and ${D}_{B}=1\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{2}\text{\hspace{0.17em}}{\mathrm{s}}^{1}$. We have launched ${10}^{8}$ photons and randomly divided the detected ones into 26 groups (to match the number of groups used in the experiments) to obtain a collection of autocorrelation functions for evaluating their variability with respect to the sampling noise and used for random sampling the IRF of our experimental system. The high number of photons was chosen to obtain good convergence of the simulated reflectance curves. We note that in order to simulate realistic experimental noise in the autocorrelation functions, a proper TDDCS noise model is needed that takes the correlator structure into account in a similar manner done for CWDCS.^{15} We have considered two broad temporal gates: an early gate, starting well before the rise of the DTOF to where it falls below 2% of its peak value (i.e., from $1$ to 0.9 ns, where we set $t=0$ at the peak of the IRF), and a late gate, from that point to where the DTOF is smaller than 0.1% of its peak value (i.e., from 0.9 to 5.9 ns). After computing the autocorrelation functions using Eq. (1), we have performed a fit for ${D}_{\mathrm{B}}$ using Eqs. (3) and (7) for an uncorrected and IRFcorrected fit, respectively. The fit minimizes the mean squared error between the simulated ${g}_{1}$ and the model, using the MATLAB (MathWorks) function “fminsearch,” which exploits the Nelder–Mead simplex method. We have chosen a fitting region from $\tau =1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{s}$ to the point were ${g}_{1}=0.7$, for the late gate, and ${g}_{1}=0.5$, for the early gate. The starting and ending points of the fit are values typically used in TDDCS experiments, to reject afterpulsing effects and regions that are effected by the probe contact, superficial heterogeneities and high noise, respectively.^{6}^{,}^{9} The different ending limit for the late gate was chosen to reject contaminations from small path length photons, which are characterized by slower decay rates and thus decorrelate at a larger $\tau $. Figure 2 shows a comparison between the theoretical autocorrelations, using the IRFcorrected and uncorrected models, and the simulated ones in the presence of the IRF, averaged over the 26 groups. In addition, a comparison between the ${D}_{B}$ retrieved with an uncorrected and IRFcorrected fit of the simulated autocorrelations is shown, also for the case of ungated acquisition (from $1$ to 5.9 ns). From Fig. 2, it is possible to see that the uncorrected model does not describe the behavior of the gated autocorrelation functions correctly, in particular in the late gate case. On the other hand, the corrected model better matches the simulated autocorrelations for both early and late gates. By making use of the corrected model for a fit of the simulated autocorrelations, the estimated ${D}_{B}$ becomes very close to the expected value ($1\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{2}\text{\hspace{0.17em}}{\mathrm{s}}^{1}$) for both the gates and also ungated acquisition. In the late gate case, the one more affected by the IRF, the error is reduced from $34\pm 1\%$ to $3\pm 2\%$. 4.1.Simulations for Different Optical Properties and SourceDetector SeparationsTo further validate the model, the simulation has been repeated for different values of optical properties and source–detector separations. We simulated six media with different combinations of absorption and reduced scattering coefficients, as specified in Table 1. In each simulation, ${10}^{8}$ photons were launched and recollected at three different SD separations: $\rho =1.2$, 2, and 3 cm. The simulated value of ${D}_{B}$ was kept to $1\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{2}\text{\hspace{0.17em}}{\mathrm{s}}^{1}$ in all cases. After the simulations were launched, the gated autocorrelations (early and late) were computed using the same temporal limits of the previous section. Figure 3 shows the fitted values of the BFI for the two temporal gates and ungated acquisition, using the uncorrected and IRF model. Table 1Simulated values of reduced scattering coefficient, absorptions coefficient, and Brownian diffusion coefficient for the six simulated media.
Using the uncorrected model, the discrepancy between early and late gate BFI is high, especially for small SD separations and high absorption coefficients, probably due to the smallest width of the path length distribution. The dependence of the error on the reduced scattering coefficient is lower. Conversely, by using the IRFcorrected model, it is possible to recover a good agreement between the gated BFIs for all the considered values of optical properties and SD separation, with an error always smaller than 3%. 5.ExperimentsWe have then performed a homogenous phantom experiment, using the liquid mixture with optical properties specified in Sec. 3.3 and an SD separation $\rho =1.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$. We have considered five temporal gates, in which the first one is centered at the peak of the IRF (i.e., $t=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ns}$). Each gate had a width of 1 ns and overlapped 500 ps to the previous one. Using the measured arrival times and ToFs, for each gate, we have computed 20 autocorrelations ${g}_{1}(\tau )$ with 1s integration time for 20 s and averaged them, repeating this 26 times, for a total of 520 s of acquisition. The group duration was chosen to increase the SNR of the autocorrelation curves while maintaining a sufficient number of groups to obtain reliable statistics from the fitted BFI. For all the autocorrelations, the fitting interval was chosen from $\tau ={10}^{6}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ to the point where ${g}_{1}=0.5$. Figure 4 shows the measured autocorrelation functions averaged over all the temporal windows and the corresponding fits, uncorrected and IRFcorrected, for an early and a late gate, centered at $t=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ns}$ and $t=1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ns}$, respectively. We also show a comparison between the retrieved BFI using the uncorrected and IRFcorrected models, together with errorbars indicating their standard deviations along the 20 s blocks, and the average ungated value retrieved fitting the corresponding autocorrelations with the CW solution of the correlation diffusion equation.^{1} As it is possible to see from Fig. 4, the uncorrected and IRFcorrected models give similar fit quality in the considered regions. The obtained value of the $\beta $ parameter decreases gradually from 0.35 for the first gate to 0.22 for the last. The $\beta $ parameter changes depending on the gate position, as already reported in Refs. 5 and 6, due to the different degree of interference arising from the gating of the DTOF and the finite laser coherence length.^{14} The tail in the measured late gate ${g}_{1}(\tau )$ at large $\tau $ may be due to larger noise induced by the lower number of photons collected at later ToFs, further study is required to clarify this point. However, the impact of the tails on the results is negligible since they were excluded by the fitting area. After the IRF correction, the discrepancy between the ${D}_{B}$ retrieved from the gated autocorrelations and the ungated ones is strongly reduced. For the later gates ($t\ge 1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ns}$), which show a larger distortion, the error is reduced from 37% to 3% on average. Since the medium is homogeneous, we do not expect any difference between earlier and later gates. Therefore, the IRF correction largely compensates the discrepancy in the earlylate result caused by the nonideal acquisition system. 6.Discussion and ConclusionIn this work, we have proposed a model to describe a realistic TDDCS experiment, characterized by a finite temporal resolution, based on measurable or known quantities, such as the IRF and the time/path length gate limits. We have introduced a function, called EGF, that quantifies the path length selection capabilities of the experimental system. In this framework, the distortions in the autocorrelation functions are explained as contaminations from path lengths in the proximity of the gate. In particular, the use of this model for IRF correction permitted to reduce the error in the retrieved BFI for the late, broad gate from 34% to 3% when MC simulations are considered. Moreover, in a homogeneous phantom experiment, the use of the proposed model reduced the BFI difference between gated and ungated autocorrelations from 37% to 3%. Regarding the limits of validity of our method, it is worth noting that the model is applicable to any gate width since it takes into account the time interval considered for the computation of the autocorrelation function. A further improvement for our method could be to consider the finite coherence length effects as in Ref. 14 by adapting the theory described in Ref. 16. While the coherence effects are not taken care of, the results in Fig. 4 and in Sec. 7 (Fig. 5) demonstrate that as long as the laser is sufficiently coherent, as indicated by relatively high $\beta $ value, the model is able to use wide gates and account for the IRF effects. In Sec. 7, we show that the model is able to correctly retrieve the dynamic properties down to a pulsewidth of 200 ps. Further modeling of coherence effects should improve this but it is beyond the scope of this work due to the complexity of accurate measurements of the coherence length of this versatile, tunable laser. We note that the method was checked for homogeneous media, but it may be generalized to heterogeneous samples, for instance, to a layered medium. To conclude, the method described in this work, which can also be applied to in vivo data, paves the way to an accurate depthresolved BF measurement. This is a critical step for many biomedical applications, where the discrimination of deeper BF from superficial BF is necessary. In addition, this work can provide useful guidance for hardware design since it permits one to understand the role of different parameters, such as detector diffusion tail, on the overall system performance. 7.Appendix7.1.Derivation of the Expression for the Nonideal Electric Field Autocorrelation FunctionIn this appendix, we report the calculations to derive Eq. (6) from Eq. (5). Starting from Eq. (5) and substituting the expressions of $\tilde{f}(s)$ and ${\tilde{g}}_{1,s}(\tau ,s)$, we obtain Eq. (12)$${g}_{1}(\tau )={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\tilde{f}(s){\tilde{g}}_{1,s}(\tau ,s)\mathrm{d}s\phantom{\rule{0ex}{0ex}}={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\frac{\tilde{R}(s)}{{\int}_{\infty}^{+\infty}\tilde{R}(s)\mathrm{d}s}\frac{{\int}_{\infty}^{+\infty}{g}_{1,s}(\tau ,s{s}^{\prime})IRF({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}}{{\int}_{\infty}^{+\infty}IRF({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}}\mathrm{d}s\phantom{\rule{0ex}{0ex}}={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}\frac{\tilde{R}(s)}{{\int}_{\infty}^{+\infty}R(s)\mathrm{d}s}\frac{{\int}_{\infty}^{+\infty}{g}_{1,s}(\tau ,s{s}^{\prime})IRF({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}}{\tilde{R}(s)}\mathrm{d}s.$$In the above passages, we have used the relations $\tilde{R}(s)={\int}_{\infty}^{+\infty}\mathrm{IRF}({s}^{\prime})R(s{s}^{\prime})\mathrm{d}{s}^{\prime}$ and ${\int}_{\infty}^{+\infty}\tilde{R}(s)\mathrm{d}s={\int}_{\infty}^{+\infty}R(s)\mathrm{d}s$; this last is valid if ${\int}_{\infty}^{+\infty}\mathrm{IRF}(s)\mathrm{d}s=1$. Then, Eq. (12) can be rewritten as follows: Eq. (13)$${g}_{1}(\tau )={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}{\int}_{\infty}^{+\infty}{g}_{1,s}(\tau ,s{s}^{\prime})IRF({s}^{\prime})\frac{R(s{s}^{\prime})}{{\int}_{\infty}^{+\infty}R(s)\mathrm{d}s}\mathrm{d}{s}^{\prime}\text{\hspace{0.17em}}\mathrm{d}s\phantom{\rule{0ex}{0ex}}={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}{\int}_{\infty}^{+\infty}{g}_{1,s}(\tau ,s{s}^{\prime})IRF({s}^{\prime})f(s{s}^{\prime})\mathrm{d}{s}^{\prime}\text{\hspace{0.17em}}\mathrm{d}s\phantom{\rule{0ex}{0ex}}={\int}_{{s}_{0}}^{{s}_{0}+\mathrm{\Delta}s}{\int}_{0}^{\infty}{g}_{1,s}(\tau ,{s}^{\prime})IRF(s{s}^{\prime})f({s}^{\prime})\mathrm{d}{s}^{\prime}\text{\hspace{0.17em}}\mathrm{d}s.$$In the above passages, we have used the relations $f(s{s}^{\prime})=\frac{R(s{s}^{\prime})}{{\int}_{\infty}^{+\infty}R(s{s}^{\prime})\mathrm{d}s}=\frac{R(s{s}^{\prime})}{{\int}_{\infty}^{+\infty}R(s)\mathrm{d}s}$ and the commutative property of convolution. Note that in the last step, we started the integration from ${s}^{\prime}=0$ since the ideal path length distribution $f({s}^{\prime})$ is equal to zero for negative path lengths. 7.2.Coherence Length Effects on Amplitude and DecayRate of the AutocorrelationsTo study the effects of limited coherence length in the autocorrelation functions, we have performed a phantom experiment decreasing gradually the width of the laser pulse. Through increasing the radiofrequency power of the acoustooptic modulator, we obtained five different values of IRF FWHM, from $\sim 500\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ps}$ down to 150 ps. We have prepared a calibrated liquid phantom with the same receipt used in Sec. 5 and used an SD separation $\rho =1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$, measuring 30 ungated autocorrelations with 10s integration time each, for a total acquisition time of 300 s. Figure 5 shows the resulting intensity and, using the Siegert relation, electric field autocorrelation functions. Also, we show the dependence of the measured $\beta $ and the estimated BFI, using the CW solution of the correlation diffusion equation, on the IRF FWHM. Decreasing the pulse width degrades the value of $\beta $ as expected by the lowering of the coherence length. However, down to IRF FWHM of the order of 200 ps, the shape of the ${g}_{1}$ curve and the BFI estimation are not significantly distorted compared to the initial configuration (maximum IRF FWHM). Thus, for pulse widths above that threshold, neglecting the coherence length term should be a good approximation, especially for gate widths smaller than the width of the path length distribution, where $\beta $ is larger compared to the ungated case. AcknowledgmentsThe research was funded by Fundacio CELLEX Barcelona; Ministerio de Economia y Competitividad/FEDER (PHOTODEMENTIA, DPI201564358C21R), Institutes de Salud Carlos III/FEDER (MEDPHOTAGE, DTS16/00087), the “Severo Ochoa” Program for Centers of Excellence in R&D (SEV20150522), the Obra Social “la Caixa” Foundation (LlumMedBcn), Institució CERCA, AGAURGeneralitat (2017 SGR 1380), LASERLABEUROPE IV, the European Union’s Horizon program (LUCA, 688303H2020ICT2015), and the EU Marie Curie Initial Training Network “OILTEBIA” (317526). ReferencesT. Durduran et al.,
“Diffuse optics for tissue monitoring and tomography,”
Rep. Prog. Phys., 73 076701
(2010). https://doi.org/10.1088/00344885/73/7/076701 RPPHAG 00344885 Google Scholar
D. J. Pine et al.,
“Diffusing wave spectroscopy,”
Phys. Rev. Lett., 60 1134
–1137
(1988). https://doi.org/10.1103/PhysRevLett.60.1134 PRLTAO 00319007 Google Scholar
F. Martelli et al.,
“There’s plenty of light at the bottom: statistics of photon penetration depth in random media,”
Sci. Rep., 6 srep27057
(2016). https://doi.org/10.1038/srep27057 SRCEC3 20452322 Google Scholar
A. Pifferi et al.,
“New frontiers in timedomain diffuse optics, a review,”
J. Biomed. Opt., 21 091310
(2016). https://doi.org/10.1117/1.JBO.21.9.091310 JBOPFO 10833668 Google Scholar
J. Sutin et al.,
“Timedomain diffuse correlation spectroscopy,”
Optica, 3
(9), 1006
–1013
(2016). https://doi.org/10.1364/OPTICA.3.001006 Google Scholar
M. Pagliazzi et al.,
“Time domain diffuse correlation spectroscopy with a high coherence pulsed source: in vivo and phantom results,”
Biomed. Opt. Express, 8
(11), 5311
–5325
(2017). https://doi.org/10.1364/BOE.8.005311 BOEICL 21567085 Google Scholar
D. Irwin et al.,
“Influences of tissue absorption and scattering on diffuse correlation spectroscopy blood flow measurements,”
Biomed. Opt. Express, 2 1969
–1985
(2011). https://doi.org/10.1364/BOE.2.001969 BOEICL 21567085 Google Scholar
A. G. Yodh, P. D. Kaplan and D. J. Pine,
“Pulsed diffusingwave spectroscopy: high resolution through nonlinear optical gating,”
Phys. Rev. B, 42 4744
–4747
(1990). https://doi.org/10.1103/PhysRevB.42.4744 Google Scholar
M. Pagliazzi et al.,
“In vivo time gated diffuse correlation spectroscopy at Quasinull sourcedetector separation,”
Opt. Lett., 43
(11), 2450
–2453
(2018). https://doi.org/10.1364/OL.43.002450 OPLEDP 01469592 Google Scholar
J. Li et al.,
“Analytical models for timedomain diffuse correlation spectroscopy for multilayer and heterogeneous turbid media,”
Biomed. Opt. Express, 8 5518
–5532
(2017). https://doi.org/10.1364/BOE.8.005518 BOEICL 21567085 Google Scholar
J. M. Tualle, E. Tinet and S. Avrillier,
“A new and easy way to perform timeresolved measurements of the light scattered by a turbid medium,”
Opt. Commun., 189
(4–6), 211
–220
(2001). https://doi.org/10.1016/S00304018(01)010458 OPCOB8 00304018 Google Scholar
J. M. Tualle et al.,
“Timeresolved diffusing wave spectroscopy beyond 300 transport mean free paths,”
J. Opt. Soc. Am. A, 23 1452
–1457
(2006). https://doi.org/10.1364/JOSAA.23.001452 JOAOD6 07403232 Google Scholar
D. Borycki, O. Kholiqov and V. J. Srinivasan,
“Reflectancemode interferometric nearinfrared spectroscopy quantifies brain absorption, scattering, and blood flow index in vivo,”
Opt. Lett., 42 591
–594
(2017). https://doi.org/10.1364/OL.42.000591 OPLEDP 01469592 Google Scholar
X. Cheng et al.,
“Time domain diffuse correlation spectroscopy: modeling the effects of laser coherence length and instrument response function,”
Opt. Lett., 43
(12), 2756
–2759
(2018). https://doi.org/10.1364/OL.43.002756 OPLEDP 01469592 Google Scholar
C. Zhou,
“Invivo optical imaging and spectroscopy of cerebral hemodynamic,”
University of Pennsylvania,
(2007). Google Scholar
T. Bellini et al.,
“Effects of finite laser coherence in quasielastic multiple scattering,”
Phys. Rev. A, 44 5215
–5223
(1991). https://doi.org/10.1103/PhysRevA.44.5215 Google Scholar
G. Maret and P. E. Wolf,
“Multiple light scattering from disordered media. The effect of Brownian motion of scatterers,”
Z. Phys. B Condens. Matter, 65
(3), 409
–413
(1987). https://doi.org/10.1007/BF01303762 Google Scholar
A. Dalla Mora et al.,
“Fastgated singlephoton avalanche diode for wide dynamic range near infrared spectroscopy,”
J. Sel. Top. Quantum Electron., 16 1023
–1030
(2010). https://doi.org/10.1109/JSTQE.2009.2035823 Google Scholar
D. Contini et al.,
“Effects of timegated detection in diffuse optical imaging at short sourcedetector separation,”
J. Phys. D Appl. Phys., 48 045401
(2015). https://doi.org/10.1088/00223727/48/4/045401 Google Scholar
L. Wang, S. L. Jacques and L. Zheng,
“MCML—Monte Carlo modeling of light transport in multilayered tissues,”
Comput. Methods Programs Biomed., 47
(2), 131
–146
(1995). https://doi.org/10.1016/01692607(95)01640F CMPBEK 01692607 Google Scholar
L. Cortese et al.,
“Liquid phantoms for nearinfrared and diffuse correlation spectroscopies with tunable optical and dynamic properties,”
Biomed. Opt. Express, 9
(5), 2068
–2080
(2018). https://doi.org/10.1364/BOE.9.002068 BOEICL 21567085 Google Scholar
M. Diop and K. S. Lawrence,
“Deconvolution method for recovering the photon timeofflight distribution from timeresolved measurements,”
Opt. Lett., 37 2358
–2360
(2012). https://doi.org/10.1364/OL.37.002358 OPLEDP 01469592 Google Scholar
D. Waithe et al.,
“FoCuSpoint: software for STED fluorescence correlation and timegated single photon counting,”
Bioinformatics, 32
(6), 958
–960
(2016). https://doi.org/10.1093/bioinformatics/btv687 BOINFP 13674803 Google Scholar
BiographyLorenzo Colombo received his MS degree in engineering physics from Politecnico di Milano, Italy, in 2017. Currently, he is enrolled in the PhD program in physics at Politecnico di Milano. His research interests focus on photon migration in diffusive media, and the development of diffuse optical systems and methods for noninvasive monitoring of biological tissues. Marco Pagliazzi is a PhD candidate in the medical optics group at ICFO, Spain. He graduated in physics at the University of Pisa, Italy. The focus of his research is developing diffuse optical instrumentation, algorithms, and protocols, and applying them to the study of human physiology. Sanathana Konugolu Venkata Sekar obtained his PhD (cum laude) from Politecnico di Milano, Milan, Italy, in 2017. Currently, he is a postdoc at Biophotonics@Tyndall National Institute (IPIC), Cork, Ireland. His current area of research and interest includes upconverting nanoparticles for deep tissue imaging and optogenetics applications, GAs in scattering media absorption spectroscopy (GASMAS) for infant lung vitals monitoring, timedomain diffuse optics, diffuse correlation spectroscopy, Raman spectroscopy, imaging and nonimaging probes development, and entrepreneurial activities in biophotonics. Davide Contini received his MS degree in electronic engineering and his PhD in physics from the Politecnico di Milano, Italy, in 2004 and 2007, respectively. He is an associate professor in the Department of Physics, Politecnico di Milano. He is the author of more than 100 papers in international peerreviewed journals and conference proceedings. His research activity is focused on timeresolved spectroscopy of highly diffusive media for applications in biology and medicine. Alberto Dalla Mora graduated summa cum laude in electronics engineering from the Politecnico di Milano, Italy, in 2006 and received his PhD summa cum laude in information and communication technology from the same university in 2010. Since 2011, he has been an assistant professor in the Department of Physics, Politecnico di Milano. He has authored more than 30 papers in international peerreviewed journals. Currently, his research interests include timeresolved diffuse spectroscopy instrumentation and applications for biomedical diagnosis. Lorenzo Spinelli graduated in physics from the University of Milan, Italy, in 1994. In 1999, he received his PhD in physics from the University of Milan. Since November 2001, he has been a researcher at the National Research Council in service at the Institute of Photonics and Nanotechnologies—Milan section. His research field has moved toward the study of photon propagation in turbid media, in order to develop optical systems for medical diagnostics. Alessandro Torricelli received his MS degree in electronic engineering from Politecnico di Milano, Italy, in 1994 and his PhD in physics from Politecnico di Torino in 1999. He is a full professor in the Department of Physics, Politecnico di Milano. He is the author of more than 100 papers in international peerreviewed journals. His current research interests include photon migration in diffusive media, fNIRS, and noninvasive diffuse spectroscopy with timedomain systems. Turgut Durduran was trained at the University of Pennsylvania. In 2009, he moved to ICFO—The Institute of Photonic Sciences, Spain, where he leads the medical optics group. His research interests revolve around the use of diffuse light to noninvasively probe tissue function. The group develops new technologies and algorithms and routinely translates them to preclinical, clinical, and industrial applications. Antonio Pifferi received his MS degree in nuclear engineering from Politecnico di Milano, Italy, in 1991 and his PhD in physics from Politecnico di Torino in 1995. He is a full professor in the Department of Physics, Politecnico di Milano. His research is directed toward the development of laser techniques and instrumentation for diagnosis and the study of light propagation in diffusive media, with applications to optical biopsy, optical mammography, and functional brain imaging. 