Simultaneous measurement of localized diffusion and flow using optical coherence tomography

We report on the simultaneous and localized measurements of the diffusion coefficient and flow velocity based on the normalized autocorrelation function using optical coherence tomography (OCT). Our results on a flowing solution of polystyrene spheres show that the flow velocity and the diffusion coefficient can be reliably estimated in a regime determined by the sample diffusivity, the local flow velocity, and the Gaussian beam waist. We experimentally show that a smaller beam waist results in an improvement of the velocity sensitivity at cost of the precision and accuracy of the estimation of the diffusion coefficient. Further, we show that the decay of the OCT autocorrelation due to flow depends only on the Gaussian beam waist irrespective of the sample position with respect to the focus position.


Introduction
Quantification of diffusive and translational dynamics of particles is interesting for the study of fundamental fluid dynamic processes such as, shear dependent diffusion [1], and for a number of applications such as e.g., aerosols [2], particle sorting [3], intracellular transport [4], biofilm growth [5], and blood circulation [6].
In principle two techniques are available to quantify flow and diffusion simultaneously. First, particle tracking based techniques such as particle tracking velocimetry (PTV), use a sequence of (microscopy) images to track the motion of individual particles which are then used to calculate the particle displacement to find the relevant dynamic parameters. Second, dynamic light scattering (DLS) based technique detect the fluctuations of single scattered light by an ensemble of particles without resolving the individual particles. By fitting an appropriate model to the autocorrelation function of the scattered light, the ensemble averaged particle dynamics are measured. However, in the case of PTV the focus has to be mechanically scanned along the propagation direction of the light to produce useful imaging ranges and in DLS the path-length distribution of the scattered light is not known, providing only volumetric averaged information about the ensemble particle dynamics.
Localized particle dynamics can be probed using optical coherence tomography (OCT). Using a combination of coherent and confocal gating OCT measures the exact path-length distribution of the scattered light up to a few millimeters deep into a sample and with high spatio-temporal resolution [7]. Functional extensions based on the combination of OCT and DLS have been developed to quantify particle diffusion [8,9,10,11], flow velocity [12,13], and flow and diffusion [14,15].
In this manuscript we report on simultaneous and localized measurements of the diffusion coefficient and flow velocity of a colloidal suspension. We quantify the fluctuations in the measured backscattered intensity by fitting a model to the normalized autocorrelation function. Previously, we have shown that the autocorrelation function can be used to quantify the flow velocity for arbitrary flow angles in the presence of diffusion [15]. In the case of shear flow the velocity gradient over the coherent detection gate is a source of additional decorrelation of the OCT signal. With the aim of studying the coupling of the flow and diffusion process in the measurement model, we constrain the measurement to be perpendicular to the flow direction. Although usually, both processes are assumed to be independent sources of decorrelation, we show that the flow velocity and the diffusion coefficient can only be reliably estimated in a regime determined by the sample diffusivity, the local flow velocity, and the Gaussian beam waist.

Optical coherence tomography system
The experiments are performed with a home built fiber-based swept-source OCT system. A schematic of the experimental set-up is shown in Fig. 1. The description of the system has been reported elsewhere [15] but is repeated here for convenience. The system operates at a center wavelength of 1312 nm with a bandwidth of 92 nm and a sweep frequency of 50 kHz (Axsun Technologies). The average output power is 20.9 mW and the duty cycle is 59.4%. Data is sampled (ATS9350, AlazarTech) with an interferometrically derived external clock signal at equidistant wavenumber intervals. To ensure phase stability each sweep is triggered by the signal of a fiber Bragg grating centered at 1266 nm (OE Land) [16]. The interferometric signal is detected with a 150 MHz balanced photodetector (PDB450C, Thorlabs) and a 80 MHz low-pass filter (VLF-80+, Mini-Circuits). The trigger signal is detected with a 125 MHz pho-  Figure 1: Schematic of the experimental swept-source OCT set-up. PD: photodetector, FBG: fiber Bragg grating, PC: polarization controllers, C: collimating lens, F: focusing lens, and FC: flow channel. Gravity is in the z-direction. Adapted from Ref. [15].
todetector (1811, New Focus). To assess the influence of the focusing optics on the measured parameters, the sample and reference arms' optics are composed of a collimating lens (PAF-X-18-C, Thorlabs) and three different achromatic doublet focusing lenses (AC254-030-C, AC254-040-C, AC254-100-C, Thorlabs). Unless otherwise stated in the text the 40 mm focal length lens (AC254-040-C) is used. The power ratio of the sample and reference arms is 90/10. The axial resolution is 8.1 ± 0.3 μm in air measured with a mirror reflector.

Flow system
Flow is generated by a perfusion pump (Perfusor fm, Braun) and directed through a cylindrical glass channel with a measured inner diameter of 1097 ± 25 μm or a glass channel with an inner diameter of 50 ± 5 μm (VitroCom). The flowing suspension consists of 1 vol. % polystyrene spheres (PPs-0.2, Kisker) dissolved in distilled water. The measured mean sphere diameter is 152 ± 15 nm (qNano, np400, Izon). The measured refractive index of the medium at 1312 nm is n = 1.33. The flow channel is placed perpendicular to the propagation direction of the imaging beam by minimizing the Doppler shift. The experimental flow conditions throughout the manuscript are well described by Poiseuille flow with a maximum Reynolds number of 9.

Data analysis
The path-length resolved diffusion coefficient and flow velocity are measured by fitting a model of the normalized autocorrelation function of the OCT signal for every path-length independently using a modified version of the model presented in Ref. [15]: where the exponential term describes the longitudinal diffusive dynamics and the Gaussian term describes the transverse directional dynamics, with z representing the optical path-length (OPL) and |τ | the time. The additional factor of 2 in the exponents accounts for the use of the magnitude of the OCT signal in the analysis [9]. D(z) is the path-length resolved diffusion coefficient given by the Stokes-Einstein equation D(z) = k B T K /6πηr, with k B Boltzmann's constant, T K the absolute temperature, η the viscosity, and r the hydrodynamic particle radius. The absolute value of the scattering vector is q = 4πn sin (α/2)/λ, with n the refractive index of the medium, λ the wavelength in vacuum, and α the scattering angle. Further, v(z) is the path-length resolved transverse flow velocity and w is the beam waist ((1/e) radius of the field). The normalization was taken with respect of g(z, 0). For our OCT set-up, the spread of q over the bandwidth is small, therefore we set q = q c at the center wavelength and α = 180 • . For the derivation of Eq. 1 we have assumed no number fluctuations, independence of particle concentration, and single scattering [15].

Model fitting
Processing of the data is performed as follows: raw interferometric data consisting of 1088 data points is Fourier transformed to calculate the complex-valued OCT signal. For every path-length we calculate the autocovariance of the magnitude of the OCT signal over 5000 time-adjacent acquisitions. This process is repeated 20 times and averaged. The transverse velocity and diffusion coefficient are determined in the time domain by fitting Eq. (1) to the normalized autocovariance of the data. In this fit, v(z) and D(z) are the only free running parameters. In all plots throughout the manuscript, all values of v and D are mean values over 5 measurements and the error bars are the corresponding standard deviations. The error on the fitted parameters is expressed as the coefficient of variation, which is defined as the ratio of the standard deviation to the mean. Throughout the manuscript the results of the model fitting are shown in the frequency domain. Figure 2 shows results for the measurement of the diffusion coefficient for a solution of polystyrene spheres in the absence of flow. Figure 2(a) shows a typical plot of the normalized OCT magnitude for the flow channel. The focus was placed near the center of the channel to reduce the influence of multiple scattering at the longer path-lengths [15]. agreement with the expected diffusion coefficient of 2.8 ± 0.3 μm 2 /s. It should be noted that boundary effects on the diffusion are to be expected, however the range of distances at which these effects occurs is small compared to the coherence length of the laser source [10]. The power spectral density at pathlengths of 516 μm and of 994 μm is calculated and plotted in Fig. 2(c) and (d), respectively. As can be observed, the data (circles) and the model (line) are in good agreement. The small deviation form the model at the end of the frequency range is attributed to noise.

Flow and diffusion, one parameter fit
We start by showing that the decay of the Gaussian term in the normalized OCT autocorrelation function in Eq. (1) is characterized by the beam waist and not by the local beam radius irrespective of sample position from the focus. Figure 3 (a) shows plots of the power spectral density for an experiment where a 50 μm diameter channel was translated in the longitudinal direction away from the position of the focus. The markers show the data corresponding to different distances from the focus, but with the same flow velocity. The data is compared to a model that depends only on the beam waist at the focus and to a model that is based on a depth dependent beam radius. As can be clearly seen, all power spectra overlap and it can be observed that the decay of the power spectral density does not depend on the value of the beam radius at the corresponding distance from the focus, but depends only on the beam waist.
The deviations from the model for high frequencies is attributed to the influence of noise which increases with increasing distance from the focus. Figure 3 (b) shows the expected beam radius for a focused Gaussian beam for a lens with a focal length of 40 mm and the values for the beam waist calculated by fitting the data shown in (a) and by letting w be the only fitting parameter in Eq.
(1). Clearly, the fitted beam waist values are independent of the position of the sample with respect to the focus. Figure 4 shows results for the measurement of the path-length resolved flow velocity with the diffusion coefficient fixed to the value measured in the no-flow case (cf. Sec. 3.1). Since the diffusion coefficient does not depend on optical path-length, v(z) is the only fitting parameter in Eq. (1). Figure 4(a) shows the path-length resolved flow velocity. The gray parabola shows the reference velocity calculated based on the flow rate set by the perfusion pump and the diameter of the channel. As can be observed, the computed flow velocity is in good agreement with the reference values. Figure 4(b) shows the coefficient of variation of the fitted flow velocity. For optical path-lengths close to the walls of the channel an increase in the coefficient of variation is observed. Figure 4(c) and (d) show plots of the power spectral density for the path-lengths shown by the arrows in Fig. 4(a). The circles show the data and the squares show the data corresponding to the no-flow diffusion experiment. For the path-length close to the wall of the channel, the decay of the power spectral density for the flow experiment is well described by a Lorentzian decay as was the case for the no-flow experiment. However, for the path-length corresponding to the center of the channel, this is no longer the case: the contribution of the Gaussian (flow) term is clearly observed in the first half of the frequency range. For all  Fig. 2(b)). (d) Similar to (c), but for a path-length corresponding to the center of the flow channel. For the sake of visualization, only every other data point has been plotted.
path-lengths, a excellent agreement is observed between the data and the model. Figure 5 shows results for the dependency of the velocity uncertainty from Gaussian beam waist. Figure 5(a) shows a plot of the beam radius for three focusing lenses measured with a knife edge [17]. The squares correspond to a lens with a focal length of 30 mm, the circles for a lens with a focal length of 40 mm, and the triangles for a lens with a focal length of 100 mm. The dashed lines show the corresponding Gaussian beam models. The measured waist values are 8.5 ± 0.5 μm, 10.8 ± 0.5 μm, and 26.3 ± 1 μm, respectively. Figure 5(b) shows the flow velocity's coefficient of variation measured in three different experiments with the three lenses. For each experiment the maximum reference flow velocity is 10 mm/s. As can be observed, the coefficient of variation monotonically decreases for increasing flow velocity for all three lenses. For a particular flow velocity, the choice of a smaller waist results in a lower value for the coefficient of variation compared to the values measured with a larger waist.   Figure 6 shows plots of the power spectral density for an optical path-length near the edge of the channel ( Fig. 6 (a)) and for an optical path-length corresponding to the center of the channel ( Fig. 6 (b)). The circles show the data and the blue line shows the fitted model. The red line shows only the contribution of the flow term and the green line shows only the contribution of the diffusion term. It can be seen that at optical path-lengths with relatively low flow velocity values the decay of the power spectral density is well described by a Lorentzian (diffusion) decay. However, for relatively larger flow velocities, the decay of the power spectral density at low frequencies is well described by a Gaussian (flow) decay and at high frequencies it converges towards a Lorentzian decay. Figure 7 shows results for the simultaneous measurement of the diffusion coefficient and the flow velocity based on data similar to the one shown in the previous figure. Figure 7(a) shows the path-length resolved flow velocities for a set of three varying reference velocities. The markers show the measured velocities and the lines show the reference velocities. As can be seen, the measured velocity values are well in agreement with the expected reference values. Figure  7(b) shows the corresponding coefficient of variation for the fitted velocities. For optical path-lengths close to the walls of the channel an increase in the coefficient of variation is observed similar to the case for a fit with fixed diffusion (cf. Fig. 4). Figure 7(c) shows the path-length resolved diffusion coefficient corresponding to the flow velocities shown in Fig. 7(a). The markers show the measured diffusion coefficient and the line shows the diffusion coefficient measured in the no-flow case which have has used here as a reference. As can be observed, the diffusion coefficient measured at relatively low flow velocities is in good agreement with the reference values. For the relatively large flow velocities at the center of the flow channel, a path-length dependency of the measured diffusion coefficient is observed where the diffusion coefficient at the center of the channel is larger than the diffusion coefficient measured near the wall of the channel. As the maximum flow velocity through the center of the channel increases, the overestimation of the diffusion coefficient increases as well. Figure  7(d) shows the coefficient of variation for the fitted diffusion coefficient. For optical path-lengths close to the center of the channel an increase in the coefficient of variation is observed. Figure 8 shows results for the simultaneous measurement of the diffusion coefficient and the flow velocity at high flow rates measured with two lenses

Discussion
At hand of the theoretical model presented in Ref. [15] we have studied the simultaneous measurement of the local diffusion coefficient and flow velocity of a colloidal suspension.
We have constrained the experimental and fitting conditions to study the influence of the diffusion coefficient, the flow velocity, and the beam waist on the measured OCT autocorrelation function and the estimated dynamical parameters. Our experimental results show that, although the diffusive and translation dynamics enter Eq. (1) independently, a reliable estimation of both quantities is challenging.
Interestingly, interpreting the Gaussian (flow) term as being a transit time effect of the moving scatterers through the illumination profile would result in a dependence on the local beam radius. However, we have shown that the decay of the normalized OCT autocorrelation function is characterized by the beam waist irrespective of the position of the sample with respect of the focus. This is in excellent agreement with the theoretical prediction given in Ref. [18], where this effect is attributed to the phase curvature of a focused Gaussian beam. More specifically, away from the focal plane the effect of beam radius increase and curvature increase compensate each other. As a result, the decay due to flow does not depend on the local Gaussian beam radius.
Our experimental results show that the uncertainty with which either the diffusion coefficient and the flow velocity can be quantified is determined by a competition between the single exponential (diffusion) and Gaussian (flow) decay terms in the autocorrelation function. In the particular case of relatively low flow velocities the decay of the autocorrelation function is dominated by the diffusive dynamics and is well described by a single exponential decay. In this case, the contribution of the Gaussian (flow) term to the total decay is small and therefore the estimation of the flow velocity suffers of a large uncertainty. In the case of relatively high flow velocities the decay of the autocorrelation function is dominated by the translational dynamics and is well described by a Gaussian decay. This results in a large uncertainty of the diffusion coefficient since the contribution of the single exponential decay to the total decay is relatively small.
The choice of the beam waist determines in which regime either dynamical process dominates the decay of the autocorrelation function. We identify two regimes for the choice of the beam waist. For applications where the focus lies in the estimation of the diffusion coefficient, it is more suitable to choose for a larger beam waist such that the decay due to diffusion dominates over a large part of the velocity range of interest. However, for applications where the focus lies in quantifying small flow velocities, a smaller beam waist is more suitable such that the decay due to flow dominates. As an illustrative example, we show in Fig. 9 the decay time constants τ D = (2Dq 2 ) −1 (diffusion) and τ v = ( √ 2v/w) −1 (flow) and the corresponding uncertainties based on the data shown in Fig. 8. We can clearly see that as long as the decay due to diffusion dominates over the decay due to flow, the uncertainty of the diffusion coefficient remains close to 5%. However, as the decay due to flow starts dominating, the uncertainty increases. A similar behavior is observed for the uncertainty of the flow velocity, where an uncertainty of 5% is measured when the decay due to flow determines the total decay. Finally, by choosing either a large or a small beam waist we can shift the regime where either, the decay due to flow or due to diffusion, determines the total decay of the autocorrelation function.
In the experiments performed for this work, we have constrained the flow to be perpendicular to the propagation direction of the imaging beam. By doing this, we have been able to study the influence of flow and diffusion on the amplitude terms of the normalized autocorrelation function. In essence, a generalization to arbitrary flow directions is possible. As we have shown previously, a longitudinal flow velocity component can be accurately measured by fitting an additional phase term to the normalized autocorrelation function [15]. For moderate shear rates the phase and the amplitude terms in the autocorrelation function remain uncoupled and thus, the results presented in this work can be generalized directly to the more complex case of arbitrary flow direction.
In principle, a generalization of our results to arbitrary values of the beam waist is possible but should be performed with care. For small beam waists a strong reduction in the depth of focus would become a limitation for practical applications requiring long imaging ranges. Also, as a consequence of tight focusing a stronger dependency on the scattering angle is expected which has been shown to result in a more complex functional decay for the diffusive dynamics [19]. Furthermore, a smaller beam waist will limit the maximum flow velocity that can be measured. In this case, the temporal resolution of the autocorrelation function should be sufficiently small to sample the flow decay τ v = ( √ 2v/w) −1 sufficiently. For large beam waists the influence of low signalto-noise ratios [13] and multiple scattering [20] should also be considered. A thorough description of these effects on the estimation uncertainty of the dynamical parameters is subject of future research.

Conclusion
We presented measurements of the local diffusion coefficient and flow velocity based on the normalized autocorrelation function measured by optical coherence tomography. Based on our model, we have obtained accurate results by fitting the model to the measured data with only the physically relevant parameters as fitting parameters. Our results on a flowing solution of polystyrene spheres show that the flow velocity and the diffusion coefficient are reliably estimated. The regime where this is possible is determined by the diffusion coefficient of the sample, the local flow velocity, and the Gaussian beam waist produced by the focusing optics. We have experimentally shown that a smaller beam waist results in an improved estimation of the flow velocity and that a larger beam waist results in an improved estimation of the diffusion coefficient. Finally, we showed that the decay of the autocorrelation due to flow depends only on the Gaussian beam waist irrespective of the sample position with respect to the focus position.