1.

Introduction

To characterize the dynamics of the flowing suspension of microscopic particles, a number of optical techniques can be used including optical coherence tomography (OCT), low coherence interferometry, and dynamic light scattering. The statistics of radiation scattered from particles undergoing purely Brownian motion in the absence of flow are well-known and established.¹ The physical reason responsible for the optical field fluctuations is the random Doppler shifts produced by the Brownian motion of the scattering particles. Similarly, the effects of purely translational flowing liquid motion (ignoring the Brownian component) on the characteristics of the scattered coherent radiation have also been examined,² yielding a dynamic speckle pattern. This translational speckle motion in the image plane of an optical system allows for the possibility of flow velocity measurement via spatial filtering. However, the quantitative statistical model of scattered radiation has not yet been developed.

Accounting for both flow and Brownian motion effects is more challenging, but is important for many OCT applications including Doppler optical coherence tomography,³ OCT imaging techniques based on dynamic light scattering,⁴ and speckle variance optical coherence tomography.^5,6 In particular, it was suggested^7,8 that the longitudinal and transverse components of flow velocity can be measured simultaneously by using the Doppler shift and the spectral width of the Doppler component, respectively. The model of spectral width used in Refs. 7 and 8 was based on empirical considerations and did not include some parameters (experimental and Brownian-motion related) which have considerable impact on the spectral width. Yet, untangling the flow and Brownian motion contributions would be desirable, in that both contain important biomedical information. For example, the statistical characteristics of red blood cells undergoing Brownian motion can potentially be used for the noninvasive measurement of blood viscosity, and hence blood glucose levels.⁹ In a living organism, the blood is seldom stagnant and it is, therefore, important to take into account the blood flow movement for in vivo measurement of blood viscosity and scattering dynamics.

Statistical models of radiation scattered by flowing Brownian particles in OCT conditions have been suggested recently.^10,11 In these papers, the solutions corresponding to the free space optical geometry of scattering have been obtained; the effects of the optics on the backscattered radiation have not been taken into account. However, most real-life OCT systems use an imaging geometry where the scattering volume is imaged onto the detector or fiber end face. Unlike random Doppler shifts, which govern the spectral broadening due to Brownian motion and do not depend on the specific OCT geometry of illumination and detection, the translational flow motion of particles gives rise to a dynamic speckle pattern. It is well known that the statistical properties of dynamic speckles are dependent on the geometrical parameters of illumination and observation.¹² Therefore, the statistical model of dynamic light scattering in the presence of flow should include the parameters of the optical system used.

In the present paper, we develop a simple analytical model for the temporal correlation function and power spectrum of optical radiation scattered by flowing Brownian particles in the image plane of a typical OCT setup comprised of two sample arm lenses. The model takes into account single scattering events only, which is justified for small penetration depths¹³ or larger penetration depths within dilute suspensions. Further, unlike the well-established model of dynamic light scattering¹ which is typically applied to relatively large scattering volumes (of the order of a cubic millimeters),¹ in the context of OCT we are considering the effects of microscopic scattering volumes on the statistical properties of scattered radiation.

2.

Correlation Function of Radiation Scattered by Uniformly Flowing Brownian Particles

Consider a collection of Brownian particles suspended in the flowing liquid under the experimental conditions of OCT. The geometry of the OCT setup is shown in Fig. 1. Here, the scattering volume is defined by the Gaussian profile of the illuminating laser beam ($x$, $y$, or lateral extent) and its coherence length ($z$ or depth extent). This volume is also called the point spread function of OCT or OCT voxel. The following equation defines the electric field of the illuminating Gaussian beam incident on the sample:

Fig. 1

(a) Optical coherence tomography (OCT) setup and (b) illumination and observation geometry of scattered radiation in OCT sample arm.

The presence of the coherence factor $\exp [-(z/{(l_{\mathrm{c}}/2)}^2]$ is due to the coherent gating effect in OCT. The Gaussian shape of it is an approximation; in fact, it may be considered as the first element of Hermite functions¹⁴ expansion of the source laser spectral density.

It is assumed that the amplitude of the reflected field is proportional to the incident field; furthermore, only single scattering events are taken into account. In the case of discrete scattering centers, the backscattered optical field can be found by using the first Born approximation for the optical field in the turbid medium¹⁵ as

Assuming stationarity, the spatio-temporal correlation function of the optical field in the image plane is given by the following ensemble average:

Putting Eqs. (1) and (3) into Eq. (2) and then substituting the result into Eq. (5) yields an expression to be summed over two indices. Since each particle moves independently (in other words, the particles do not interact with each other), the averages, which contain different indices, yield zero. It means that the summation over two indices can be replaced by a summation over only one index.^10,11 The ensemble average of each term of the sum is performed by convolving it with the probability density

Assuming that the number of particles in the scattering volume is large, the summation over discrete particles can be replaced by integration over the initial positions of the particles ${\mathit{r}}_0$. The evaluation of the Gaussian integral thus obtained is straightforward yielding the following result for the normalized temporal correlation function:

The equivalent Gaussian beam radius $w_{\mathrm{e}}$ is defined as

Strictly speaking, Eq. (10) is valid only in the image plane of a two-lens optical system since the Green’s function used [Eq. (3)] is valid only here. Therefore, the dependence of $w_{\mathrm{e}}$ on the displacement from the image plane [$w=w(z)$, $\rho =\rho (z)$] from Eq. (10) can be viewed only as a contribution of Gaussian beam radius and wavefront curvature change into the dependence of $w_{\mathrm{e}}(z)$. The effects of defocusing can be modeled through the complex ABCD matrix formalism¹⁶ which does not yield a simple analytical solution. The corresponding numerical calculation that takes into account all three factors predicts that the equivalent Gaussian beam radius is independent of the position of scattering volume relative to the Gaussian beam waist $w_{\mathrm{e}}(z)=w_{\mathrm{e}}(z=0)=\mathrm{const}$.

Quite similarly, we are obtaining the correlation function in the case of free space by using the following expression for Green’s function:¹²

The correlation function in the free space has the same appearance as Eq. (8), the only difference being that the equivalent Gaussian beam radius is here

At the Gaussian beam waist $\rho =\infty$ and the term containing $kw/[2(z-Z)]$ can be neglected as it is much smaller than the $1/w$ term. In typical OCT conditions with $w=10\mu \mathrm{m}$, $\lambda =1.3\mu \mathrm{m}$, and $z-Z=100\mathrm{mm}$, we get $kw/[2(z-Z)]=0.24{\mathrm{mm}}^{-1}$, $1/w=100{\mathrm{mm}}^{-1}$. Although the term containing $kw/[2(z-Z)]$ is small, it cannot be entirely neglected in the range of displacements along the $z$-axis where the wavefront curvature $\rho$ takes its minimum $\rho \sim 2z_{\mathrm{F}}$ and where $z_{\mathrm{F}}$ is the Gaussian beam Rayleigh range.

So, in the free space we get $w_{\mathrm{e}}\approx w$ with the scattering volume at the Gaussian beam waist; in the image plane of a two-lens optical system we get $w_{\mathrm{e}}\approx w/\sqrt{2}$.

Equation (8) has a structure similar to one previously reported in Ref. 11, but with important differences. In the typical OCT backscattering imaging geometry, the dependence on $F_{\mathrm{s}}$ and $q$ does drop out and the only important difference between our two-lens imaging geometry and the free space situation is the different value of the equivalent Gaussian beam radius. However, there is some value in seeing the full equation without any approximations to note the relative importance of dependencies on various system parameters; further, for a transmission-type OCT system one has to use Eq. (10) with parameters $F_{\mathrm{s}}$ and $q$ included, as the fiber and/or the optics used to collect the scattered radiation are typically different from the emitting fiber or optics shaping the incident laser beam. There may be other scenarios (e.g., cross-talk considerations in multichannel OCT) that may also benefit from the insights afforded by the full equation of Eq. (10).

To interpret the result in Eq. (8), we note that the factors $s_{\mathrm{w}}^{-1}$ and $s_{\mathrm{c}1}^{-1}$ in front of the exponential terms are caused by speckle fluctuations due to Brownian motion. The first exponential term is the contribution of random Doppler shifts due to Brownian motion. The second exponential term represents a Doppler shift due solely to flow motion. The other Brownian motion factor $s_{\mathrm{c}2}$ clarifies that there is a frequency modulation in the Doppler component. The third exponential term is the contribution of speckle fluctuations due to translational flow.

Equation (8) thus shows that in the general case, the Brownian motion and translational flow motion contributions do not appear as independent factors (i.e., the second and third exponential “flow” terms have admixtures of Brownian motion contributions). This means that these two movements cannot be untangled or evaluated independently, at least in the context of the full expression of Eq. (8). However, we can make several OCT-relevant approximations (below) and examine the resultant regimes where these two contributions do indeed separate and can be evaluated independently.

The impact of the factors given in Eq. (9) can be estimated from the following considerations. The correlation time due to Brownian motion is smaller than the correlation time due to flow in the low velocity range. Therefore, here, we have to estimate the impact of factors Eq. (9) over the time $\sim \tau =s_{\mathrm{c}1}{(2k^{2}D)}^{-1}\approx {(2k^{2}D)}^{-1}$. For the mean square displacement, we get

The terms with the spatial scale $l_{\mathrm{c}}/2$ in Eq. (9) may have importance only in the low velocity range when using a very broadband source with a submicron coherence length.¹⁷ The term containing the laser beam radius (typically $w_{\mathrm{e}}\sim 5\mu \mathrm{m}$ or so) may be smaller on the order of a few microns in OCT systems with dynamic focus.¹³ Thus, in most practical OCT situations, we have

By using these simplifying approximations which are valid in most OCT regimes, Eq. (8) reduces to following equation for the correlation function of scattered radiation:

As a simple check of the derived equation, note that in the absence of flow, Eq. (11) reduces to a known dynamic light scattering result for the correlation function of the optical field backscattered by Brownian particles:¹

The correlation function of the photodetector current fluctuations is obtained by taking the real part of Eq. (11), since the photocurrent generated by the OCT photodetector is proportional to the real part of $\langle E_{\mathrm{s}}{E}_{\mathrm{r}}^{*}\rangle$ where $E_{\mathrm{r}}$ is a reference arm optical field:

In Eq. (12), ${\tau }_{\mathrm{b}}$ is a correlation time due to Brownian motion,

Note that the range of validity of Eq. (12) is not limited strictly to the optical system shown in Fig. 1, where the distance $L$ between collimating and sample lenses is equal to the sum of the focal distances $L=F_{\mathrm{c}}+F_{\mathrm{s}}$. The geometry shown in Fig. 1 was chosen because it yields a simple analytical solution for the correlation function of scattered radiation. The approach based on an ABCD matrix enables the derivation of the correlation function for an optical system of any complexity in the paraxial approximation.¹⁶ The corresponding calculation shows that within the Rayleigh range of the collimated beam $Z_{\mathrm{F}}=\pi q^2/\lambda \gg L$, the dependence of correlation time ${\tau }_{\mathrm{t}}$ on the distance between lenses $L$ is negligibly small. In fact, in this case, the radius of the collimated beam is equal to the radius of the limiting (Gaussian) diaphragm for the backscattered radiation. Thus, Eq. (12) can be used for any imaging optical system comprising a collimating and sample lens provided that $Z_{\mathrm{F}}\gg L$.

Recapping, we have derived the analytical expression for the temporal correlation function of radiation scattered by flowing Brownian particles, taking into account the specifics of a standard OCT imaging system, including the effect of the small OCT scattering volume. In general, the contributions of random Doppler shifts due to Brownian motion and speckle due to translational flow cannot be untangled (i.e., presented as independent factors in the resultant correlation function [Eq. (8)]). Indeed, this full expression would correctly describe the temporal correlation function of a high-resolution OCT system (submicron optical resolution). On the other hand, in more common OCT systems with moderate (5 to $25\mu \mathrm{m}$) spatial resolution, the contributions of Brownian motion and particle translational flow movement can be presented as independent factors in the joint correlation function [Eq. (12)] via realistic approximations.

3.

Power Spectrum of Scattered Radiation

The power spectrum of the scattered radiation can be found via a Fourier transform of the correlation function in Eq. (12). Since Eq. (12) consists of a product of a negative exponential, a Gaussian, and a cosine, the resultant power spectral shape is given by the convolution of Fourier transforms of the negative exponential and Gaussian terms centered at the Doppler frequency (arising from the Fourier transform of the cosine term). The Fourier transform of a negative exponential is a Lorentzian, while that of a Gaussian function is also a Gaussian. Thus, the resultant spectral shape of the flow-caused Doppler peak [cosine term in Eq. (12)] is given by the convolution of the Lorentzian and Gaussian functions:

The convolution of the Lorentzian and Gaussian yields a so-called Voigt function, a well-known mathematical description commonly arising in the studies of atomic and molecular spectra.¹⁸ To aid computation, one can use the relationship between the Voigt function and the complementary error function. Following Ref. 19, the convolution in Eq. (15) can be expressed as

The physical interpretation of the Voigt spectrum shape is quite different in atomic or molecular spectroscopy versus the power spectrum of scattered radiation derived here. In the former, the Lorentzian spectral broadening is due to collisions, and Gaussian broadening is caused by Doppler shifts due to thermal motion of molecules. Here, the nature of spectral shape contributions is quite different—the Lorentzian spectral shape is caused by random Doppler shifts of the Brownian motion, and the Gaussian component is caused by speckle fluctuations due to the translational flow particle movement.

The spectral width of the Doppler peak is a convenient tool for spectrum characterization, thus its quantification would be desirable. Unfortunately, since the Voigt function cannot be expressed through elementary analytical functions, one has to use various approximations to get a simple empirical dependence of the Voigt spectral width as a function of the spectral widths of a Gaussian and a Lorentzian. The following empirical equation²⁰ has been proposed to obtain the full spectral width at half maximum (${\mathrm{FWHM}}_{\mathrm{D}}$) of a Voigt function

Equations (13), (14), (19), and (20) show that apart from the transverse velocity component $v_{\mathrm{t}}=v\sin (\theta )$, the longitudinal velocity component, $v\cos (\theta )$ and the Brownian motion component contribute to the full spectral width of the Doppler peak. These additional contributions deserve careful consideration when trying to extract the transverse velocity from the experimentally measured Doppler peak width, something that may have been overlooked in earlier publications.^7,8

4.

Results and Discussion

Exploring the predictions of the developed formalism through the derived spectral width of the Doppler component in Eq. (19), we now consider several examples of OCT M-mode measurements of spherical Brownian particles suspended in a flowing stream of water. In the discussions which follow, we assume spherical $1\text{-}\mu \mathrm{m}$ diameter Brownian particles in water at 22°C interrogated with an OCT system (central wavelength $\lambda =1.31\mu \mathrm{m}$, illuminating Gaussian beam waist radius $w=10\mu \mathrm{m}$). The value of water viscosity was calculated from an empirical equation for viscosity dependence on temperature: $\eta (22°\mathrm{C})=0.955\times {10}^{-3}\mathrm{Pa}·\mathrm{s}$.²¹

Figure 2 shows the flow velocity at which the contributions to the spectral width of Brownian motion and speckle fluctuations (due to translation/flow) are equal, i.e., ${\mathrm{FWHM}}_{\mathrm{b}}={\mathrm{FWHM}}_{\mathrm{t}}$. Equations (13), (14), (19), and (20) have been used in the calculations. Several OCT scattering volumes, as determined by relative magnitudes of lateral (beam focusing) and longitudinal (coherence length) resolutions have been examined. In the range of angles and velocities below each curve, the Brownian motion dominates the statistics of scattering; conversely, in the range of $(v,\theta )$ above the curve, flow-induced speckle fluctuations prevail. Figure 2 suggests practical choices of the relevant experimental regimes: for example, for optimal measurements of Brownian motion parameters (e.g., diffusivity) in the presence of flow, one needs to operate in the $(v,\theta )$ parameter space below the curve. In the regions of $(v,\theta )$ above the curve, the contribution of flow-induced speckle fluctuations will severely affect the accuracy of diffusivity measurements. Note that the shape of the curves with the maximum or minimum at $\theta =90\mathrm{deg}$ is caused by the angle-dependent relative contributions of transverse and longitudinal velocity components to the dynamics of the speckle pattern. For a (spherically symmetric) OCT scattering volume with $l_{\mathrm{c}}=2w_{\mathrm{e}}$ (curve 2), these relative contributions are independent of the Doppler angle.

Fig. 2

The flow velocity curves at which the contributions to the spectral width of translational and Brownian motions are equal, as a function of angle between the flow velocity vector and the OCT axis. Curve 1: $l_{\mathrm{c}}=4w_{\mathrm{e}}$, curve 2: $l_{\mathrm{c}}=2w_{\mathrm{e}}$, and curve 3: $l_{\mathrm{c}}=w_{\mathrm{e}}$.

Figure 3 shows the FWHM of the Doppler peak as a function of flow velocity. The angle between the flow velocity vector and the OCT axis for this simulation was kept fixed at $\theta =45\mathrm{deg}$, and Eqs. (13), (14), (19), and (20) were used. Note that the range of velocities shown in Fig. 3 does not start at zero but rather at $v=15\mu \mathrm{m}/\mathrm{s}$. This is because a certain minimum velocity is needed for the spectral shape to exhibit both high and low frequency wings, and thus for the concept of FWHM to be meaningful. At lower velocities, the low frequency wing begins to disappear.

Fig. 3

The spectral width at half maximum of the Doppler component as a function of flow velocity for $\theta =45\mathrm{deg}$. Curve 1: full width at half maximum, curve 2: contribution of translational flow motion, and curve 3: contribution of Brownian motion.

Figure 3 demonstrates that for these simulation parameters (notably with $\theta =45\mathrm{deg}$), the presence of flow with velocities less than $\sim 100\mu \mathrm{m}/\mathrm{s}$ does not significantly modify the Brownian motion nature of the scattered radiation—the impact of speckle fluctuations is small. In other words, the Brownian particle movement contribution to the spectrum of the Doppler component is defining the dynamics of the scattered radiation. As expected, for increasingly fast flows, the contribution of the translational component begins to dwarf the Brownian motion effects.

In light of recent reports that derive the transverse flow velocity component from the broadening of the Doppler peak,^7,8,10 consider the explicit dependence of the absolute flow velocity at which the transverse flow contribution to the spectral width is equal to the joint contribution of the Brownian and longitudinal flow motions, i.e.,

Figure 4 shows the system behavior for three different regimes determined by OCT focusing/coherence length properties (in analogy with Fig. 2). The area of parameters $(v,\theta )$ where the transverse flow contribution prevails is only above each curve. In the regions below, to the left and to the right of each curve, the joint contribution of Brownian motion and longitudinal flow motions dominate and dwarf the transverse flow component contribution. This means that the reliable measurements of the flow velocity transverse component are possible only in a limited range of Doppler angles near $\theta =90\mathrm{deg}$, with the exact details determined by the shape of the OCT scattering volume.

Fig. 4

The absolute flow velocity curves at which the contribution to the spectral width of transverse flow equals the joint contribution of Brownian motion and longitudinal flow, as a function of angle between the flow velocity vector and the OCT axis. Curve 1: $l_{\mathrm{c}}=4w_{\mathrm{e}}$, curve 2: $l_{\mathrm{c}}=2w_{\mathrm{e}}$, and curve 3: $l_{\mathrm{c}}=w_{\mathrm{e}}$.

Let us now compare the developed formalism with the theoretical and experimental results previously published.^10,11 As already noted in Sec. 1, Refs. 10 and 11 present an analysis of backscattered radiation based on the solution of the scattering problem in free space. The integral equation (12) in Ref. 10, when evaluated in the absence of a flow velocity gradient, yields a simple analytical equation that closely resembles Eq. (11) in the present paper: the equivalent Gaussian beam radius is equal to $w/\sqrt{2}$, where the factor $\sqrt{2}$ is introduced empirically in Ref. 10 to account for the fiber coupling efficiency. Equation (10) of the present paper yields the same result for the equivalent Gaussian beam radius. However, here, the physical reason for the appearance of the factor $\sqrt{2}$ is the different character of the speckle motion in free space and the image plane. As follows from the analysis given in Ref. 12, the speckle movement character is predominantly “boiling” in the plane of the sample lens for the parameters of illumination and observation used in Ref. 10. At the same time, in a two lens optical imaging system (as used experimentally in Ref. 10), a translational component in speckle motion appears, which makes an additional contribution to the speckle dynamics and increases the spectral width. Mathematically, it is expressed as a smaller equivalent Gaussian beam radius: $w_{\mathrm{e}}=w/\sqrt{2}$.

It is unlikely that the coupling efficiency of the fiber alone can modify the temporal characteristics of the dynamic speckle pattern since it characterizes the optical fiber in static conditions. As per discussion leading to Eqs. (8)–(10), the temporal analysis of the backscattered radiation in the image plane of the specific optical system is needed to obtain a valid result for the temporal correlation function and the spectrum of scattered radiation.

The comparison of the flow velocities obtained from mass discharge and spectral measurements when the scattering volume is at the center of the Gaussian beam waist has shown a reasonable agreement.¹⁰ This means that the theoretical predictions given in the present paper are also experimentally validated.

In Ref. 11, the equivalent Gaussian beam radius is not introduced, so that $w_{\mathrm{e}}=w$, which applies to the geometry of free space. All the other parameters of the temporal correlation function in Eq. (11) of the present paper are the same as in the correspondent result obtained in Ref. 11—apart from the wavefront curvature, which is not taken into account in Ref. 11.

When comparing the results of Ref. 11 with our model, it is important to note that their Gaussian beam radius is not defined explicitly; instead, the authors are using “the inverse of $1/e$ width of OCT transverse resolution” denoted as $h_{\mathrm{t}}$. The transverse part of their point spread function takes the form $\exp [-2h_{\mathrm{t}}^2(x^2+y^2)]=\exp [-2(x^2+y^2)/w^2]$. This means that the definitions of $w$ in Ref. 11 and in our paper differ by $\surd 2$, since the transverse part of PSF in our paper has the form $\exp [-(x^2+y^2)/w^2]$ as follows from Eq. (1). Therefore, although the results for the correlation time of scattered radiation in Ref. 11 look identical to our model of a two-lens imaging system, they actually differ by $\surd 2$; the model in Ref. 11 does yield identical results with ours evaluated at the Gaussian beam waist in free space.

Another important issue is the presence, importance, and minimization of the so-called “pedestal” component in the spectrum of scattered radiation situated near zero frequency. There are a number of physical reasons which contribute to this component: self-beating²² of the scattered radiation, fluctuations in the particle density, instrumental $1/f$ noise, etc. In the absence of flow, the pedestal component is a significant source of noise that adds to the spectral density of the optical field fluctuations caused by Doppler shifts due to Brownian motion, and cannot be separated from it. However, under conditions of flow, the pedestal and Doppler components are separated. This allows the estimation of the pedestal in the photocurrent spectrum and its effective removal. Paradoxically then, while significantly fast flows can certainly obscure the Brownian motion contributions, slow to moderate flows can actually improve the accuracy of the Brownian motion fluid properties determination by shifting the peak with flow and Brownian motion contributions to higher frequencies away from the pedestal noise peak. For a set of experimental parameters (interrogation angle, scattering volume, etc.) one can thus potentially use the flow properties in a well-designed M-mode OCT measurement to accurately measure Brownian motion dynamics. Experimental validation of some of these ideas will be pursued in a separate publication.

5.

Conclusions

We have developed a simple analytical model for the correlation function and power spectrum of optical radiation scattered by flowing Brownian particles for an imaging geometry of a common OCT setup, and have examined various approximate regimes relevant to a typical OCT measurement. The model includes the parameters of the optical system used to shape the Gaussian beam in the sample and collect the backscattered radiation on the fiber end face. It is shown that the spectral shape of the power spectrum, centered at the flow-determined Doppler frequency, is given by a Voigt function, a mathematical descriptor arising from atomic and molecular spectral shapes. The developed formalism allows one to quantify the relative contributions of the random (Brownian) and translation (flow) components of the scattered radiation signal. Derived dependencies on the Doppler angle and flow velocity suggest experimental operating regimes where Brownian motion can be quantified even in the presence of flow. Further, we have re-examined the possibility of transverse flow velocity component determination based on the measurement of the spectral width of the Doppler shifted peak. It is shown that with a typical OCT geometry, these measurements are possible only in a limited range of Doppler angles near $\theta =90\mathrm{deg}$. Finally, it is shown that the introduced empirically equivalent Gaussian beam radius in Ref. 10 has the same value at the waist of the Gaussian beam as predicted by our analysis of scattered radiation in the image plane of a two lens optical system.