Statistical properties of dynamic speckles from flowing Brownian scatterers in the vicinity of the image plane in optical coherence tomography

Ivan Popov; Andrew Weatherbee; I. Alex Vitkin

doi:10.1364/BOE.8.002004

1. Introduction

Optical techniques based on the study of radiation dynamics from moving or deforming biotissue have received considerable attention in recent years. In particular, optical coherence tomography (OCT) allows for simultaneous localized measurements of the translational flow velocity and diffusivity of particles suspended in a flowing liquid [1,2]. Further, dynamic statistics of the optical radiation scattered in biological tissue contain important biomedical information, which can be utilized for blood velocity measurements [3] and blood microvascular network imaging [4–6]. The scattered radiation also carries information about diffusion properties of the scattering particles, which may potentially yield blood viscosity (likely dependent on the shear rate since blood is a non-Newtonian fluid) and thus a link to blood glucose levels [7]. Since blood in vivo is rarely stagnant, one has to develop the methodology to extract diffusion characteristics from the flowing suspension of Brownian particles.

Until recently, the models describing the radiation scattered from flowing Brownian particles were based on the free space geometry, i.e. they did not take into account the impact of the imaging system [1,8,9]. For example, Weiss et al [1] fit the experimentally obtained spectra to a free space model with an additional empirically-introduced factor. This is problematic in several ways, including the issue that empirical models do not have a well-defined applicability area. Therefore, it is desirable to develop closed-form analytical models of scattering which specifically include the parameters describing the imaging setup used in the experiment.

It is well known that the coherent radiation scattered by particulate media gives rise to a speckle pattern [10] and that the speckle statistics are different in the free space vs the image plane of a given optical system [11]. Since OCT yields depth resolved measurements, it is important to get the solution of the scattering problem not only for the exact focal plane of the optical system used, but also in its vicinity, specifically within several Rayleigh ranges of it. In our recent efforts, an ABCD ray tracing matrix methodology was introduced to obtain numerical solutions of the scattering problem for a number of frequently used optical imaging geometries in the OCT sample arm [12,13]. However, it is preferable to have a closed-form analytical expression for the spatio-temporal correlation function of the scattered radiation. This will yield information not only about temporal correlation time and spatial correlation radius of the scattered radiation, but will also allow analysis of the speckle motion; in other words, one can relate the speckle velocity in the image plane to flowing particle velocity in the image plane as would be observed with incoherent illumination.

In the present paper, we develop such a mathematical model for the scattered optical field in the OCT sample arm, thus obtaining a closed-form analytical expression for the spatio-temporal correlation function. The model assumes the use of a clean imaging optical system in the OCT sample arm; however, as will be shown, the developed framework is also relatively accurate for other commonly used OCT sample arm optics. We also analyze the characteristics of the dynamic speckle motion in the fiber end face plane, and derive the expressions for the speckle radius, speckle velocity, and correlation time due to speckle movement. This analysis explains the previous observation [2] that the correlation time of the backscattered field does not depend on the axial position in the scattering volume relative to the focal plane. In addition, we analyze how the OCT speckle velocity is related to the particle image velocity under incoherent illumination conditions. For experimental validation, OCT measurements in flowing particulate scattering suspensions with varying viscosities (varying amounts of dissolved glucose) were performed, demonstrating good agreement between the measured data and model predictions.

2. Statistical model of scattered radiation in the OCT fiber end face plane

Let us consider the collection of small spherical scatterers (Brownian particles) suspended in a liquid and flowing through a glass capillary. It is assumed that the flow is laminar; apart from the bulk translational flow motion, the particles are experiencing random chaotic displacements due to Brownian motion (translational diffusion).

A two-lens optical system in the OCT sample arm shapes the Gaussian beam emerging from the single-mode fiber end, so that its waist (focal plane of the optical system) is positioned in the flowing liquid (Fig. 1). The backscattered radiation is collected by the same fiber, combined with the reference beam, and directed to the OCT detector. The optical arrangement shown in Fig. 1(b) is known as a ‘clean optical system’.

Fig. 1 (a) OCT set-up. (b) Illumination and observation geometry of scattered radiation in an OCT sample arm. F_c and F_s are the focal distances of the collimating and sample lenses, respectively; z₁ is the axial offset between the scattering volume and the beam waist. The collimating and sample lenses and the Gaussian diaphragm (Gaussian variation of transparency in the radial direction) form what is called the ‘clean optical system’. Note that the diaphragm need not be physically present in the set-up, as its role is played by the acceptance angle of the fiber.

Download Full Size | PDF

The electric field of a Gaussian beam can be written as

E_{i n} (r, t) = E_{i n 0} \frac{w_{0}}{w} \exp [- (\frac{1}{w^{2}} - \frac{i k}{2 ρ}) (x^{2} + y^{2}) + i k z - i ω t + i φ_{G}] \exp [- {(\frac{z - z_{1}}{l_{c} / 2})}^{2}] .

Here r = (x,y,z) is the Cartesian coordinate system with the origin at the center of the beam waist, E_in0 is the amplitude of the electric field of the incident optical wave at the centre of the laser spot at the beam waist, w is the radius of the Gaussian beam at 1/e amplitude level, w₀ is the beam waist radius, k = 2πn/λ₀ = 2π/λ, λ is the wavelength in the medium, λ₀ is the wavelength in air, n is the refractive index of the medium, ρ is the radius of wavefront curvature, ω is the laser beam angular frequency, and ϕ_G is the Gouy phase [14]. The last exponential term in Eq. (1) is due to the OCT coherence gating effect; in this term l_c is the coherence length of the illuminating laser, and z₁ is the position of the scattering volume relative to the beam waist. The familiar Gaussian form of this term is an approximation, and can be viewed as a zeroth-order term of the Hermite functions expansion [15] of the Fourier transform of the source laser spectral density.

As shown in [12]_, in the single scattering regime, the spatio-temporal correlation function C_s of radiation scattered from flowing Brownian particles (point scatterer model) can be represented as the product of the contributions of Brownian motion C_sb and translational motion due to flow C_st (provided that l_c/2, w₀ >> λ/(2π)):

C_{s} (X_{1}, Z, X_{2}, Z, τ) = 〈 E_{s} (X_{1}, Z, t_{1}) E_{s}^{*} (X_{2}, Z, t_{2}) 〉 = C_{s b} (τ) C_{s t} (X_{1}, Z, X_{2}, Z, τ)

where <…> denotes ensemble averaging, the asterisk sign * denotes complex conjugation, τ = t₁ – t₂, and X, Z are coordinates in the fiber end face plane with X = (x,y). The correlation is between the temporal points t₁ and t₂ at spatial points X₁ and X₂.

The contribution due to Brownian motion is well known [16]:

C_{s b} (τ) = \exp (- {(2 k)}^{2} D_{d} τ)

where D_d is particle diffusivity.

The contribution due to translational flow motion can be represented as [11]:

\begin{array}{l} C_{s t} (X_{1}, Z, X_{2}, Z, τ) = \\ \int E_{i n} (x, z, 0) K (x, z, X_{1}, Z) E_{i n}^{*} (x + v_{x y} τ, z + v_{z} τ) K^{*} (x + v_{x y} τ, z + v_{z} τ, X_{2}, Z) d x d y d z \end{array}

where x = (x,y); X = (X,Y); v_xy = (v_x,v_y), v_z are flow velocity components, and K(r,R) is Green’s function of the optical system used. The current formalism assumes the single scattering regime, that the number of particles in the scattering volume is large, that the conditions for a fully developed speckle pattern are met [17,18], and that the scattering system is ergodic (i.e., ensemble averaging is assumed equal to the time averaging).

In general, the Green’s function of any optical system in the paraxial approximation is given by the following equation [19]:

K (r, R) = - \frac{i k}{2 π B} \exp {\frac{- i k}{2 B} (D X^{2} - 2 x X + A x^{2}) + i k (z - Z)}

where A, B, and D are (generally complex-valued) elements of the ABCD matrix. The following expressions for the matrix elements of the optical system shown in Fig. 1(b) have been obtained in [13]:

A = - \frac{F_{c}}{F_{s}}, B = - \frac{F_{c} z_{1}}{F_{s}} - \frac{2 i F_{c} F_{s}}{k q^{2}}, C = 0, D = - \frac{n F_{s}}{F_{c}}

where q is the radius of the Gaussian diaphragm at 1/e² intensity, and F_c and F_s are shown in Fig. 1(b). The diaphragm is used in the mathematical model to allow for the fiber acceptance angle effect to be taken into account. In practice, its experimental presence is optional in that its role is played by the acceptance angle of the fiber. In the model, the diaphragm radius is equal to the radius of the collimated Gaussian beam at the 1/e² intensity level, which in turn is determined by the numerical aperture of the fiber emitting the Gaussian beam:

q \approx N A_{f} F_{c} \approx \frac{λ_{0} F_{c}}{π w_{f}} = \frac{2 F_{c}}{k_{0} w_{f}} = \frac{2 F_{s}}{k_{0} w_{0}}

where NA_f is the numerical aperture of the fiber, w_f is the mode field radius of the fiber, and k₀ is the wavenumber in air.

By putting q = 2F_s/(k₀w₀) into the expression for the matrix element B given by Eq. (6), and then substituting the result into Eq. (5) we arrive at:

K (r, R) = - \frac{i k}{2 π (z + i z_{F})} \exp {- (\frac{1}{w^{2}} + \frac{i k}{2 ρ}) [{(x + M X)}^{2} + (n - 1) {(M X)}^{2}] + i k (z - Z)}

where M = F_s/F_c is the magnification of the optical system and z_F = πw₀²/λ is the Rayleigh range of the Gaussian beam in the medium. Here we also made use of the well-known equations for the Gaussian beam radius and wavefront curvature [14]:

w = w (z) = w_{0} \sqrt{1 + {(\frac{z}{z_{F}})}^{2}}, ρ = ρ (z) = z + \frac{z_{F}^{2}}{z} .

Equation (8) gives a somewhat unexpected result: it has largely the same appearance as Eq. (1) for the incident optical field. This is because the limiting diaphragm radius q for backscattered radiation is equal to the radius of the incident collimated Gaussian beam.

By putting Eq. (1) and Eq. (8) into Eq. (4), we get the integral equation for the spatio-temporal flow correlation function C_st. Then the integration with respect to x, y and z can be performed separately. The integral over coordinates (x,y) is a Gaussian integral, which can be evaluated in the usual way. However, the integral over z is not a Gaussian integral, as both the incident electric field of the optical wave and the Green’s function of the optical system (Eqs. (1) and (8), respectively), are showing a rather complicated dependence on z. Specifically, the Gaussian beam radius w, wavefront curvature ρ, and Gouy phase ϕ_G are z-dependent. The variation scale of these quantities is the Rayleigh range z_F.

To proceed, the integration with respect to z can be performed by noting that, when compared with the coherence gating term, the laser beam radius, wavefront curvature, and Gouy phase do not change significantly on the scale of half of the coherence length l_c/2, and can thus be assumed constant in z. Indeed, with the typical Gaussian beam radius w₀ = 10 μm, wavelength λ₀ = 1.3 μm, medium refractive index n = 1.32 (water), one gets for the Rayleigh range z_F = πnw₀²/λ₀ = 420 μm; in contrast, a typical coherence length l_c for common OCT systems is in the range 5 to 20 μm.

So, the presence of the coherence gate term in Eq. (1) greatly simplifies the z-integration in Eq. (4): provided that l_c/2<<z_F, the integral with respect to z contains only linear and quadratic terms in the exponent and can be evaluated in the usual way. Indeed, the z-dependent term in Eq. (4), obtained after substitution of Eqs. (1) and (8) into Eq. (4), takes the following form:

\exp [- \frac{2}{{(l_{c} / 2)}^{2}} {(z - z_{1} + \frac{v_{z} τ}{2})}^{2}] .

This means that the point of extremum with respect to z is at

z_{m a x} = z_{1} - \frac{v_{z} τ}{2}

and thus the laser beam radius, wavefront curvature, and Gouy phase should all be evaluated at this point.

Let us estimate the magnitude of the time-dependent term in Eq. (10). We are not interested in times exceeding the combined correlation time, which includes the correlation times associated with Brownian motion as well as the transverse and longitudinal flow components. The correlation time τ_cz due to the longitudinal flow component is larger than the combined correlation time. This is because the more physical processes there are affecting the spectral width, the more it widens; conversely, the corresponding correlation time decreases. Therefore, let us estimate τ_cz to be used as an upper bound on the combined correlation time. This is of the order of time needed for a particle to cross the OCT coherence scattering volume in the z direction, i.e.τ_cz ~l_c/v_z. This means that v_zτ_cz/2 appearing in Eq. (10) is ~l_c/2. Therefore, with l_c <<z_F, the term v_zτ/2 can be safely neglected and we can set z_max ≈z₁_.

Performing the resultant integrations with respect to (x,y) and z in Eq. (4) yields the following result for the translational part of the spatio-temporal correlation function (normalized to unity at z₁ = 0):

\begin{array}{l} C_{s t} (X_{1}, Z_{f}, X_{2}, Z_{f}, τ) = \\ {(\frac{ω_{0}}{ω})}^{2} \exp (- \frac{X_{1}^{2} + X_{2}^{2}}{2 ω^{2}}) \exp {2 i k v_{z} τ - {[\frac{X_{1} - X_{2} + \frac{2}{M} (v_{x y} τ)}{r_{s p}}]}^{2} - \frac{1}{2} {(\frac{v_{z} τ}{l_{c} / 2})}^{2}} \end{array}

where

r_{s p} = \frac{2 w_{0}}{M} = 2 w_{f}

is the speckle radius. In other words, the speckle radius is twice the mode field radius of the fiber. Also in Eq. (11),

ω = \frac{w_{0}}{M \sqrt{2 n - 1}} \sqrt{1 + \frac{z_{1}^{2}}{z_{F}^{2}}}

is the radius of the backscattered radiation beam in the fiber end face plane with the position of the scattering volume shifted by z₁ from the Gaussian beam waist.

Note that apart from the intensity term containing the beam radius ω, Eq. (11) does not show any dependence on z₁, the position of the scattering volume relative to the beam waist. This is because after integration of Eq. (4), the resulting combination of the Gaussian beam radius and wavefront curvature cancels out the dependence on z₁. Therefore, in this single-scattering point scatterer model [12], the correlation properties of the spatio-temporal correlation function (the second exponential term in Eq. (11)) do not depend on the position of the scattering volume relative to the beam waist. In addition, Eq. (12) explicitly shows that the speckle radius also does not depend on the position of the scattering volume relative to the beam waist z₁.

Equation (12) also shows that speckle radius does not depend on the parameters of the optical system or the properties of the particulate scattering medium. However, similar calculations for surface scattering [20] show that for the so-called random continua medium (continuous profile of the refractive index n(r), as is the case for biological tissue), the speckle radius will depend on the parameters of the medium, and can thus potentially provide information about the interrogated tissue. The case of modelling such random continua media will be addressed in a future publication.

As an interesting aside, with X₁ = X₂ = X the exponential term in Eq. (11) containing ω is the intensity profile of the backscattered radiation beam in the fiber end face plane: $I (X, z_{1}) = {(\frac{ω_{0}}{ω (z_{1})})}^{2} \exp (- {(\frac{X}{ω (z_{1})})}^{2})$ , where ω₀ is the backscattered beam radius with z₁ = 0. Also, with X₁ = X₂ = 0, τ = 0, Eq. (11) yields the ratio of backscattered beam intensities I(z₁)/I(0) = ω₀/ω(z₁)². This result allows estimation of SNR degradation for points in OCT imaging planes offset from the beam waist (z₁ = 0).

With X₁ = X₂ = 0, i.e. in the centre of the fiber end face, Eqs. (11), (12), and (13) yield the following for the temporal correlation function of Eq. (2):

C_{s} (τ) = C_{s b} (τ) C_{s t} (τ) = \exp (2 i k v_{z} τ) \exp (- \frac{τ}{τ_{b}}) \exp [- {(\frac{τ}{τ_{t}})}^{2}]

where

τ_{b} = {[{(2 k)}^{2} D_{d}]}^{- 1}

is the correlation time due to Brownian motion (as per Eq. (3)), and

τ_{t} = {[\frac{v_{x}^{2} + v_{y}^{2}}{w_{0}^{2}} + \frac{1}{2} {(\frac{v_{z}}{l_{c} / 2})}^{2}]}^{- 1 / 2} = {v^{2} [{(\frac{\sin (θ)}{w_{0}})}^{2} + \frac{1}{2} {(\frac{\cos (θ)}{l_{c} / 2})}^{2}]}^{- 1 / 2}

is the correlation time due to dynamic speckle caused by translational flow motion.

The first two terms in Eq. (14) are well-known: the first is the Doppler shift term, and the second is the contribution due to Brownian motion. However, the third term is less familiar. It is the contribution of dynamic speckles due to the flow motion, with τ_t given by Eq. (16). The same equation appears also in Eq. (7) of Ref [13], and in Eq. (11) of Ref [12]. In Ref [12], this expression for τ_t was derived analytically at only the position of the scattering volume at the Gaussian beam waist; and in Ref [13], using numerical methods, the same expression was shown to be valid not only at the beam waist but also within several Rayleigh ranges of it. The present paper uses a purely analytical approach, and the derived expression for τ_t is valid for any displacement from the beam waist. It is also valid for many practical OCT sample arm imaging setups (provided these remain paraxial optical imaging systems), whereas the validity of Eq. (7) of Ref [13] has only been demonstrated for selected parameter values used in the numerical procedure (as per Table 2 of Ref [13].).

In Eq. (14), the second (Brownian motion) term shows negative linear time dependence in the exponent, whereas the third term (dynamic speckle fluctuations due to flow) shows negative quadratic dependence. The physical reason for the decorrelation of the Brownian contribution is the displacement of the Brownian scatterers in time via the process of diffusion (mean squared displacement), on the scale of the wavelength. The physical reason for the decorrelation of the flow-caused speckle fluctuations component is that the scatterers are passing into and out of the Gaussian profile of the illuminating laser beam; this profile is what leads to the Gaussian decorrelation term in Eq. (14) which contains the quadratic dependence on time in its negative exponent.

Let us analyze the speckle motion characteristics due to flow by examining Eq. (11). For this, we introduce the correlation coefficient defined by the following equation:

\begin{array}{l} c_{s t} (X_{1}, Z_{f}, X_{2}, Z_{f}, τ) = \\ \frac{C_{s t} (X_{1}, Z_{f}, X_{2}, Z_{f}, τ)}{\sqrt{I (X_{1}, z_{1})} \sqrt{I (X_{2}, z_{1})}} = \exp {2 i k v_{z} τ - {[\frac{X_{1} - X_{2} + \frac{2}{M} (v_{x y} τ)}{r_{s p}}]}^{2} - \frac{1}{2} {(\frac{v_{z} τ}{l_{c} / 2})}^{2}} . \end{array}

The analysis of speckle motion based on the correlation coefficient given by Eq. (17) shows that purely transverse flow (v_z = 0) causes purely translational speckle movement, even with the OCT scattering volume outside of the focal plane. Indeed, with v_z = 0 and $X_{2} - X_{1} + \frac{2}{M} v_{x y} τ = 0$ , we have c_st = 1, which means there is direct correlation between spatial points X₁ and X₂. It also means that the speckles are moving according to $X_{s p} (τ) = X_{0} - \frac{2}{M} v_{x y} τ$ , where X₀ is the speckle position at τ = 0; the speckle velocity is $v_{s p} = - \frac{2}{M} v_{x y}$ . Note the presence of the factor of 2 in this equation: the speckle velocity is 2x higher than the velocity of the scattering particles seen in the image plane under incoherent imaging conditions; similarly, speckle displacement is 2x larger than the displacement of the scattering particles as seen in the image plane. Also note that, like the speckle radius in Eq. (12), the speckle velocity does not depend on the out-of-focal plane location of the scattering volume. It follows from Eq. (17) that translational flow motion along the z-axis (v_xy = 0) produces pure speckle “boiling”. In other words, individual speckles deform, disappear, and reappear, without appreciable displacement of their positions [11]. Indeed, with v_xy = 0 the motion equation for speckles given above becomes $X_{s p} (τ) = X_{0}$ ; the speckles are not moving.

Also, with v_xy = 0 the remaining terms in Eq. (17) do not depend on z₁, meaning the speckle statistics due to “boiling” are independent of the position of the scattering volume with respect to the beam waist. We can now see all of the physical reasons for the correlation time being independent of the position of the scattering volume relative to the focal plane: the speckle radius, the speckle velocity, and the statistics of speckle “boiling” all do not depend on the position z₁ of the scattering volume relative to the focal plane. The result that the speckle velocity and displacement are 2x higher than the particle image velocity and displacement, as seen in the image plane under incoherent imaging conditions, might be important in various contexts, for example in displacement based OCT elastography [21,22].

It should be clarified that the model developed in this paper cannot be directly applied to the speckle statistics of 3D OCT images, due to the specifics of the various methods (galvo mirrors, etc.) used in practice to perform the OCT scan along the x, y, and z axes. While outside the scope of this paper, the analysis of 3D OCT image speckles should, in the future, be based on the present formalism supplemented with a separate study to account for (i) the specific 3D scanning patterns, and (ii) the fact that the statistics of scattered radiation engendered by scanning are not always equivalent to the statistics in the case of translational movement of the sample [23].

The power spectrum of scattered radiation is obtained as a Fourier transform of the temporal correlation function given by Eq. (14). As shown in [12], this spectrum can be described via the Voigt function by the following expression:

W (f) = 2 \sqrt{π} τ_{b} U (2 π τ_{b} f, \frac{τ_{b}}{τ_{t}})

where

U (x, t) = Re [\sqrt{\frac{π}{4 t}} e^{z^{2}} e r f c (z)], z = \frac{1 - i x}{2 t}

and U(x,t) is the Voigt function [24].

Recapping, by using the Green’s function approach, we have obtained the expression in Eq. (11) for the translational spatio-temporal correlation function for backscattered radiation in the fiber end face plane. This equation contains the information about correlation time and spatial speckle radius; its analysis yielded speckle motion characteristics in the fiber end face plane. We now turn to OCT experiments, with the aim to validate the developed theoretical model.

3. Experimental setup

A research grade fiber based spectral domain (SD-OCT) system was used, operating in M-mode for this study. The detailed description of the system is given in [25] and the main points are highlighted briefly here for convenience. The SD-OCT system is powered by a super-luminescent diode unit (D-1300 Superlum Ltd.) with emission bandwidth of 113 nm at FWHM, centered at 1313 nm and with 12 mW optical power. The detection end of the system is comprised of a high resolution spectrometer (P&P Optica Inc.), interfaced with a 1024 pixel, linear array IR CCD camera (SU1024LDH-1.7µm, Sensors Unlimited Inc.), with a readout rate selected to be 1123 Hz for all experiments described below.

In the sample arm, the radiation is emitted from the single mode fiber and shaped by an optical system, which produces a Gaussian beam with 11.5 μm at 1/e² intensity level radius at the waist in air. The laser beam profile was characterized using a scanning slit optical beam profiler (Thorlabs Inc.). The axial OCT resolution defined by the source laser coherence length l_c was measured by using a front-surface mirror, yielding l_c = 13 μm full width at 1/e OCT signal level in air.

A microbore glass capillary with an inner diameter of 165 μm (Accu-glass Inc.) was placed perpendicularly to the OCT optical axis in the sample arm. The choice of perpendicular geometry was to avoid the effects due to velocity gradients in the capillary on the measurement results. In the non-perpendicular geometry, the spectrum will be shifted, and the resultant velocity gradient within the OCT scattering volume will cause spectral broadening. The simple theory developed in Section 2 does not take these effects into account; future work will address this.

The scattering phantoms consisted of 215.1-nm-diameter polystyrene microspheres (Bangs Labs Inc.) suspended in water with varying amounts of dissolved D-glucose to modulate viscosity. The microsphere concentration was 0.5% by volume in all phantoms. One phantom contained spheres suspended in water, while two additional phantoms had spheres suspended in 800 mM and 1600 mM glucose solutions, respectively. The resulting scattering coefficients, calculated via Mie theory, were 0.24 cm⁻¹ (glucose-free), 0.21 cm⁻¹ (800 mM), and 0.18 cm⁻¹ (1600 mM); this variation in μ_s for a constant concentration of scattering microparticles is due to the glucose’s refractive index matching effect [26].

The flow was driven by a syringe pump actuated by a stepper motor operating in 17 nm steps (New Era Pump Systems Inc.), with an air-tight glass barrel 100uL volume syringe (Hamilton Co.). In the reported experiments, flow velocity was either zero (no-flow suspensions), or 1.94 mm/s at the capillary centre (flowing suspensions). Operating in M-mode, for each A-scan, 1024 spectrum amplitude values were obtained equally spaced in wavelength, then resampled via interpolation to produce 1024 amplitudes equally spaced in wavenumber. Taking the inverse Fourier transform of these 1024 wavenumber amplitudes produced the complex OCT signal with 512 complex values, as a function of depth into the sample. Each M-mode data set had 0.5 mega-samples of data for each depth. In this study we only analyzed the data taken from the center of the capillary. For the center of capillary depth, the power spectrum of the scattered radiation was found by taking the square modulus of the Fourier transform of the real part of the complex OCT signal.

Since the temporal correlation function, and correspondingly the power spectrum, do not depend on the displacement z₁, there was no need to accurately position the beam waist within the capillary (such as placing the beam waist at the center of the capillary, or any other position of interest).

4. Results and discussion

Commonly in speckle research, it is seen that the correlation time of the optical field fluctuations, due to the transverse velocity of the sample, depends on the speckle radius and the characteristics of speckle motion such as speckle translation and “boiling” [11]. These characteristics are difficult to predict from qualitative considerations, since they depend not only on the parameters of the optical system but also on the beam spot radius and the wavefront curvature in the scattering volume. Therefore, our derived result that the correlation coefficient is independent of displacement of the scattering volume from the beam waist (and thus also independent of the Gaussian beam radius and wavefront curvature) is noteworthy: many parameters impacting this function cancel each other, yielding independence of the speckle radius and the correlation time on the axial displacement of the scattering volume.

Similar independence of correlation time on the axial displacement was previously observed in a homodyne (non-interferometric, single-beam) experiment under perpendicular geometry [27]; however, the mathematical model in that study is based on the free-space geometry and, strictly speaking, cannot be applied to imaging set-ups. As we have shown [13], imaging and free space geometries do indeed produce different results for the correlation time of radiation scattered from flowing particles in OCT conditions. Under an OCT setup that does not have any optics in the sample arm (free space imaging conditions), the analysis of speckle motion due to transverse flow shows that there results pure “boiling” speckle motion [11]. In contrast, under the imaging geometry considered in the present paper, transverse particle flow produces pure translational speckle motion (no “boiling”). This is the physical reason which explains the difference in correlation times between the free space and imaging geometries. We also note that longitudinal particle flow produces pure boiling in both the free space and OCT imaging geometry setups.

The analytical result given by Eq. (14) echoes our earlier numerical analysis of the flow contribution to the correlation function based on ABCD ray tracing matrices given in [13]. The comparison of the mathematical model given in the present paper and the numerical model based on ABCD ray tracing matrices [13] shows nearly complete agreement for the clean imaging optical system (within computer calculation accuracy ~10⁻¹⁵). Apart from the clean imaging optical system, a number of other OCT sample arm optical geometries (one lens, two lens and multiple lenses) have been studied [13]. It was shown that for all systems with optical elements in air, the correlation time of the scattered radiation is approximately the same, with deviation from the clean imaging optical system of no more than 1%. This implies that the simple analytical results obtained in the present paper can also be used for other sample arm optical systems. In [13], we have also compared the ABCD matrix based model with the experimental results published in [2], which has confirmed the independence of the scattered radiation correlation time on displacement of the scattering volume relative to the optical system focal plane.

Now, let us experimentally validate the power spectrum functions Eqs. (18), (19), and thus the temporal correlation function Eq. (14), of the real part of the complex OCT signal arising from non-flowing and flowing Brownian particles. Since particle diffusivity and viscosity of glucose solutions (measured at the temperature of 23.8 ± 0.3 °C in our experiments) are not known from the literature, we chose to first measure them under stagnant (no-flow) conditions. Here, the power spectrum of the scattered radiation has the Lorentzian shape:

L (f) = \frac{1}{\frac{1}{{(2 π τ_{b})}^{2}} + f^{2}} .

The experiments conducted under stagnant conditions allow measurement of the correlation time τ_b by fitting Eq. (20) to the experimentally obtained power spectra, with τ_b as the fitting parameter. Both experimental data (points) and resultant fits (lines) are shown in Fig. 2. The diffusivity obtained from Eq. (15) was used to calculate the viscosity of the solution by using the Einstein-Stokes equation for spherical particle diffusivity $(k_{B} T) / 6 π η a$ with a = 107.6 ± 0.8 nm, T = 23.8 ± 0.3 °C, λ₀ = 1313 ± 0.1 nm, n = 1.32 ± 0.003. The values of solution viscosity thus found are given in Table 1, for the experimental temperature of 23.8 °C.

Fig. 2 Power spectra of backscattered radiation from stagnant Brownian particles (polystyrene microspheres) in the aqueous solutions of 0 mM, 800 mM, and 1600 mM glucose. Symbols are experimentally obtained power spectrum values; solid lines are theoretical fits of Eq. (20) with τ_b as the fitting parameter.

Download Full Size | PDF

Table 1. Correlation time, diffusivity and viscosity in stagnant suspensions of microspheres in water with varying amounts of dissolved glucose. All quoted errors were calculated following [28], and are given as the RMS error. The last italicized row shows the corresponding viscosity results measured under flow conditions (see Fig. 3 and text).

View Table

The reference viscosity of water [29] at 23.8 °C is η = 0.914 mPa·s, showing very good agreement with the experimentally obtained value of 0.90 ± 0.01 mPa·s. The small 1.5% discrepancy, outside the measurement uncertainty of ± 0.01 mPa·s, may be due to the slight diameter polydispersity of the “monodispersed” microspheres as purchased from the manufacturer.

With this data in hand, we can now proceed to the flow experiments using the same samples, to validate the theory predictions summarized in Eqs. (18) and (19). Figure 3 shows the power spectra (points) obtained by processing the real part of the complex OCT signal; the solid lines show the theoretical dependence for the spectra calculated from Eqs. (18) and (19). The correlation times τ_b and τ_t were found by using Eqs. (15) and (16) with the no-flow viscosity results given in Table 1 (3rd row); in calculating τ_t we set the experimental parameters as v = 1.94 ± 0.01mm/s, θ = π/2, w₀ = 11.5 ± 0.3 μm. We also fit Eqs. (18) and (19) to the same perpendicular flow data as shown in Fig. 3, using τ_b as the only fitting parameter (τ_t was a fixed parameter calculated using the same values of v, θ and w₀ as above); note that the solid lines in Fig. 3 are showing theoretical dependencies, not fits. From the fitted values of τ_b (relative errors of the fits were less than 1%) we calculate the experimentally measured viscosity of the three phantoms using D_d=(k_B T)/6πηa. The results are shown in the italicized bottom row of Table 1. These flow-experiment viscosities agree well and are consistent with the viscosity values from the no-flow conditions listed in the row above.

Fig. 3 Power spectra of backscattered radiation from flowing Brownian particles in solutions containing 0, 800 and 1600 mM glucose. Crosses, triangles and circles are experimentally obtained power spectrum values; solid lines are theoretical dependencies given by Eqs. (18) and (19).

Download Full Size | PDF

It is seen from Fig. 3 that there is good qualitative and quantitative agreement between model predictions (Eqs. (18) and (19)) and the experimentally obtained power spectra. This, along with the fact that the viscosity values obtained for the flow case agree well with the no-flow results (two bottom rows of Table 1), gives credence to the validity of the developed theoretical formalism. It means that the measured power spectra of scattered radiation from flowing Brownian particles provide accurate information about solution viscosity. The corresponding error analysis shows that the typical uncertainty of glucose solution viscosity measurements is ~1% in the no-flow (stagnant) experiments, and range from 4 - 10% in the perpendicular flow experiments. The main contribution to the error in the flow case is the uncertainty in the Gaussian beam radius w₀ (on the order of 3%). The results of fitting the Voigt function to the data shown in Fig. 3 has shown that the difference between the value of τ_b obtained in the no flow case and the flow case does not exceed 0.5%.

In practice, such as for bio-imaging, when both the velocity and Doppler angle are unknown, it is necessary to fit the Voigt profile to the measured spectrum with three unknown parameters: τ_b, τ_t, and the Doppler frequency. In this case, the fitted value of τ_b will still allow for the determination of the diffusion coefficient and viscosity. For the perpendicular flow data collected in the present study, with the Doppler frequency set to zero, we have performed the two-parameter Voigt fitting of τ_b and τ_t, and have found that τ_b remained constant to within ~1% across the capillary. Although it has not yet been reported in the literature, the three-parameter fit for the Doppler angle flow case will also work equally well, since it is a simple matter to include in the fitting the Doppler frequency peak position. Other recent publications have also demonstrated that it is possible to fit both τ_b and τ_t simultaneously in the perpendicular geometry set-up, without loss of accuracy in the measured τ_b. For example, the results in ref [2] shows that within a certain range of flow velocities, the accuracy of the measurement of τ_b (and the diffusion coefficient D) does not change significantly in the two parameter fit, provided that the Voigt width contribution due to translation flow (FWHM of the Gaussian) is not more than ~5x the contribution to the width by Brownian motion (FWHM of the Lorentzian), i.e. τ_t must remain greater than ~τ_b/5. In other words, when the velocity becomes too high (relative to the Brownian motion contribution), the accuracy of the measurement of the diffusion coefficient begins to suffer.

Future work will extend our methodology to examine finer viscosity variations of more direct physiological relevance (e.g., mM range for the glucose monitoring problem in the blood of diabetic patients). The additional complexities arising from the non-Newtonian properties of blood, non-spherical Brownian particles as per erythrocyte scattering, and of non-perpendicular interrogation angles, will also be examined.

5. Conclusions

We have developed a simple mathematical model for the spatio-temporal correlation function of backscattered radiation in the OCT sample arm fiber end face plane, which assumes the single scattering regime, a particulate scattering medium, and a clean optical system in the sample arm. We show that transverse flow causes purely translational speckle motion in the image plane, for any position of the OCT scattering volume relative to the focal plane (beam waist). Also, we derive the result that the OCT speckle velocity is 2x higher than the velocity of the scattering particles in the image plane under incoherent imaging conditions. Further, the speckle radius is 2x larger than the fiber mode field radius and does not depend on the parameters of either the sample arm optics or the particulate scattering medium; also, it does not depend on the position of the scattering volume relative to the focal plane. It then follows that the correlation time due to transverse translational flow also does not depend on the position of the scattering volume relative to the focal plane. Our M-mode experiments using suspensions of Brownian particles, under flowing and stagnant conditions, demonstrate good agreement between theoretical and experimentally obtained power spectra of backscattered radiation over a wide range of viscosities. Apart from applications to flow velocity and particle diffusivity quantification, the developed formalism may prove relevant elsewhere, for example in displacement based OCT elastography. In addition, this methodology for accurate viscosity determination may become useful for non-invasive glucose monitoring.

Funding

This study was supported by the Canadian Institutes of Health Research (Grant No. 126172), the Ministry of Education and Science of the Russian Federation (Grant No. 14.B25.31.0015), and the Natural Sciences and Engineering Research Council of Canada / Canadian Institutes of Health Research through the Collaborative Health Research Program (CHRP Grant No. J365581-09).

Acknowledgments

The authors thank Dr. Kostadinka Bizheva (University of Waterloo) for her help with the Fourier Domain OCT system.

References and links

1. N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Localized measurement of longitudinal and transverse flow velocities in colloidal suspensions using optical coherence tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 88(4), 042312 (2013). [CrossRef] [PubMed]

2. N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Simultaneous measurement of localized diffusion and flow using optical coherence tomography,” Opt. Express 23(3), 3448–3459 (2015). [CrossRef] [PubMed]

3. V. X. D. Yang and I. A. Vitkin, “Principles of Doppler OCT,” in Optical Coherence Tomography in Cardiovascular Research, E. Regar, T.G. van Leeuwen, P. W. Serruys eds., (Informa Healthcare, London, 2007).

4. J. Lee, W. Wu, F. Lesage, and D. A. Boas, “Multiple-capillary measurement of RBC speed, flux, and density with optical coherence tomography,” J. Cereb. Blood Flow Metab. 33(11), 1707–1710 (2013). [CrossRef] [PubMed]

5. V. Kalchenko, K. Ziv, Y. Addadi, N. Madar-Balakirski, I. Meglinski, M. Neeman, and A. Harmelin, “Combined application of dynamic light scattering imaging and fluorescence intravital microscopy in vascular biology,” Laser Phys. Lett. 7(8), 603–606 (2010). [CrossRef]

6. A. Mariampillai, B. A. Standish, E. H. Moriyama, M. Khurana, N. R. Munce, M. K. K. Leung, J. Jiang, A. Cable, B. C. Wilson, I. A. Vitkin, and V. X. D. Yang, “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett. 33(13), 1530–1532 (2008). [CrossRef] [PubMed]

7. H. Ullah, A. Mariampillai, M. Ikram, and I. A. Vitkin, “Can temporal analysis of optical coherence tomography statistics report on dextrorotatory glucose levels in blood?” Laser Phys. 21(11), 1962–1971 (2011). [CrossRef]

8. T. W. Taylor and C. M. Sorensen, “Gaussian beam effects on the photon correlation spectrum from a flowing Brownian motion system,” Appl. Opt. 25(14), 2421–2426 (1986). [CrossRef] [PubMed]

9. J. Lee, W. Wu, J. Y. Jiang, B. Zhu, and D. A. Boas, “Dynamic light scattering optical coherence tomography,” Opt. Express 20(20), 22262–22277 (2012). [CrossRef] [PubMed]

10. Y. Aizu, T. Ushizaka, and T. Asakura, “Measurements of flow velocity in a microscopic region using a transmission grating: a differential type,” Appl. Opt. 24(5), 627 (1985). [CrossRef] [PubMed]

11. T. Yoshimura, “Statistical properties of dynamic speckles,” J. Opt. Soc. Am. A 3(7), 1032–1054 (1986). [CrossRef]

12. I. Popov, A. S. Weatherbee, and I. A. Vitkin, “Dynamic light scattering arising from flowing Brownian particles: analytical model in optical coherence tomography conditions,” J. Biomed. Opt. 19(12), 127004 (2014). [CrossRef] [PubMed]

13. I. Popov and A. Vitkin, “Dynamic light scattering by flowing Brownian particles measured with optical coherence tomography: impact of the optical system,” J. Biomed. Opt. 21(1), 017002 (2016). [CrossRef] [PubMed]

14. A. Yariv, Quantum Electronics (Wiley, 1989), pp. 676.

15. Hermite function. M.V. Fedoryuk (originator), Encyclopedia of Mathematics. URL: http://www.encyclopediaofmath.org/index.php?title=Hermite_function&oldid=18370 Accessed August 24, 2016.

16. B. J. Berne and R. Pecora, Dynamic Light Scattering (Dover Publications, New York, 2000).

17. G. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser Speckle and Related Phenomena, J. C. Dainty, ed., 2nd ed., (Springer, 1984), pp. 342.

18. J. C. Dainty, “The statistics of speckle patterns,” Prog. Opt. 14, 3–46 (1976).

19. S. A. Collins Jr., “Lens-system diffraction integral written in terms of matrix optics,” J. Opt. Soc. Am. 60(9), 1168–1177 (1970). [CrossRef]

20. H. T. Yura, B. Rose, and S. G. Hanson, “Dynamic laser speckle in complex ABCD optical systems,” J. Opt. Soc. Am. A 15(5), 1160–1166 (1998). [CrossRef]

21. J. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). [CrossRef] [PubMed]

22. J. Rogowska, N. Patel, S. Plummer, and M. E. Brezinski, “Quantitative optical coherence tomographic elastography: method for assessing arterial mechanical properties,” Br. J. Radiol. 79(945), 707–711 (2006). [CrossRef] [PubMed]

23. L. M. Veselov and I. A. Popov, “Characteristics of scattered radiation under coherent-beam scanning along a rough surface,” Opt. Spectrosc. 70, 636–639 (1991).

24. N. M. Temme, “Chapter 7: Error Functions, Dawson’s and Fresnel Integrals,” in NIST Handbook of Mathematical Functions, F. W. J. Olver, ed., (Cambridge University, 2010).

25. B. Davoudi, A. Lindenmaier, B. A. Standish, G. Allo, K. Bizheva, and A. Vitkin, “Noninvasive in vivo structural and vascular imaging of human oral tissues with spectral domain optical coherence tomography,” Biomed. Opt. Express 3(5), 826–839 (2012). [CrossRef] [PubMed]

26. C. R. Hammond, CRC Handbook of Chemistry and Physics (2016).

27. T. W. Taylor and C. M. Sorensen, “Gaussian beam effects on the photon correlation spectrum from a flowing Brownian motion system,” Appl. Opt. 25(14), 2421–2426 (1986). [CrossRef] [PubMed]

28. BIPM, IEC, ISO IFCC, and IUPAP IUPAC. “OIML 1995 Guide to the Expression of Uncertainty in Measurement.” ISO, Geneva 3 (1995).

29. D. A. Berstad, B. Knapstad, M. Lamvik, P. A. Skjølsvik, K. Tørklep, and H. A. Øye, “Accurate measurements of the viscosity of water in the temperature range 19.5–25.5 C,” Physica A 151(2), 246–280 (1988). [CrossRef]

	Pure water	800 mM glucose	1600 mM glucose
τ_b, ms	2.78 ± 0.01	4.03 ± 0.01	6.15 ± 0.02
Diffusivity, μm²/s	2.26 ± 0.01	1.553 ± 0.006	1.019 ± 0.004
Viscosity, mPa•s	0.90 ± 0.01 0.914 – from literature [29]	1.30 ± 0.01	1.96 ± 0.02
Viscosity, mPa•s (flow experiment)	0.90 ± 0.04	1.30 ± 0.09	2.0 ± 0.2

Statistical properties of dynamic speckles from flowing Brownian scatterers in the vicinity of the image plane in optical coherence tomography

Abstract

1. Introduction

2. Statistical model of scattered radiation in the OCT fiber end face plane

3. Experimental setup

4. Results and discussion

5. Conclusions

Funding

Acknowledgments

References and links

Cited By

Figures (3)

Tables (1)

Equations (21)

Biomedical Optics Express