1.

Introduction

Statistical properties of coherent radiation scattered by an ensemble of flowing Brownian particles and measured by optical coherence tomography (OCT) have received considerable attention. This interest is driven by a wide range of biomedical applications for quantitative blood flow measurements,^1,2 vascular network imaging,^3–6 and intracellular motility.⁷ In addition to flow parameters, the scattered radiation also carries information about diffusion properties of scattering particles, which may potentially yield blood viscosity and thus a link to blood glucose levels.⁸ Indeed, it was demonstrated recently that the determination of flow velocity and diffusivity of Brownian particles may be possible with OCT.⁹

However, the relevant OCT models of scattered radiation developed to date (e.g., Refs. 9 and 10 and citations therein) are based on free-space solutions and do not include the parameters of optical system used to collect the scattered light. This is justified only in the absence of flow, where fluctuations of scattered radiation are caused by random Doppler shifts due to particle Brownian motion and do not depend on system optics. However, it is well known that the presence of flow gives rise to a dynamic speckle pattern,¹¹ and its statistical properties are different in the free space and at various planes of an imaging optical system.¹²

OCT signals can be processed to yield depth-resolved correlation functions and spectra of backscattered radiation. Suitable OCT mathematical models can then link parameters of Brownian motion and flow with correlation and spectral properties of scattered radiation. Thus, depth-resolved values of flow velocity vector components and particle diffusivity are obtainable. Therefore, it is important to develop the mathematical model not only in the focal plane of an OCT optical system, but also outside of it. Indeed, it is stressed in recent papers^13,14 that there is a lack of experimentally validated quantitative models for OCT measurements of flowing Brownian particles.

We have recently developed a new analytical model for correlation function and spectrum of scattered radiation to address this need¹⁵ and to supplement/extend the previous formalisms.^9,10 Since its applicability was limited to the case of the so-called 4f system (distance between collimating and sample lenses is equal to the sum of focal distances), additional study is warranted to extend its use to other useful (arbitrary) imaging geometries. This requires some novel derivations, as speckle statistics are generally dependent on first-order parameters of the optical system (e.g., lens focal distances, limiting diaphragm diameter, and distances between lenses¹²).

An ABCD matrix approach has been able to derive solution for spatiotemporal speckle statistics for any optical system in the case of surface scattering.¹⁶ However, to the best of our knowledge, the applicability of this method to volume scattering is yet to be examined. Although ABCD technique does not furnish simple insights into the impact of optical system (cf, a closed-form analytical model), it does yield a numerical solution at any plane of an optical system consisting of any number of lenses, mirrors, and media with different refractive indices (arbitrary optical properties). We thus use this approach to derive the correlation function of scattered radiation at the fiber-end face plane, for a number of practical and widely used OCT-specific sample arm geometries. The resultant model predictions are compared with existing models and with experimental data available in the literature.

2.

Mathematical Model of Coherent Radiation Scattered by the Flowing Brownian Particles

Consider an OCT optical system shown in Fig. 1. The laser beam exiting the fiber end is shaped and focused by the ABCD optical system onto the collection of suspended flowing Brownian particles and scattered back into the fiber. We have the following equation for the electric field of the incident Gaussian laser beam:

Fig. 1

OCT schematic with scattering volume shifted relative to the sample arm focal plane.

As shown in Ref. 15, if the size of scattering volume in each direction is considerably larger than $\lambda /(2\sqrt{2}\pi )$ and contains a large number of particles, the correlation function of scattered radiation is given by a product of two terms: one due to Brownian motion and one due to flow. The Brownian motion term is well known:¹⁸

The contribution of flow can be obtained by generalizing the known expression for correlation function for surface scattering to the volumetric scattering case:¹⁹

In the general case of an ABCD optical system (paraxial approximation), the Green’s function is²⁰

The integration of Eq. (3) with respect to axial variable $z$ can be performed easily if the variation scale of matrix elements $A$, $B$, and $D$ in $z$ is much larger than the source coherence length $l_c$. The estimate of this scale is given in Appendix. From there, calculations show that under typical OCT conditions, the scale of matrix element $B$ $z$-variation is $l_B\sim z_F$, where $z_F$ is the Rayleigh range of laser beam in the medium. With typical values of medium refractive index $n=1.38$, wavelength ${\lambda }_0=1.31\mu \mathrm{m}$, and Gaussian beam radius at the waist $w_0=10\mu \mathrm{m}$, we get $z_F=n\pi w^2/{\lambda }_0=330\mu \mathrm{m}$, whereas typical OCT coherence lengths are 5 to $20\mu \mathrm{m}$. The numerical estimates further show that the matrix elements $A$ and $D$ have a weak dependence on $z$, meaning that they can be treated as constants during the integration of Eq. (3) with respect to $z$ as well.

The optical field $E_{in}(\mathbf{r},t)$ in Eq. (1) contains other $z$-dependent terms, specifically $w(z)$, $\rho (z)$, and ${\phi }_G(z)$. As is known,¹⁶ the scale of variation of these terms is $\sim z_F$; thus all $z$-dependent terms in Eq. (3) [apart from coherence gate and $ikz$ terms in the incident optical field $E_{in}(\mathbf{r},t)$] are changing slowly on the scale of coherence length and can be assumed constant. The remaining Gaussian integral, which contains coherence gate term, can be evaluated in the usual way. Thus, the derivations stemming from the proposed ABCD approach are valid provided that the source coherence length is much smaller than the Rayleigh range of the beam in the sample. This condition is easily realized in most OCT systems, with the possible exceptions of highly focussed optical coherence microscopy setups.

The procedure of integrating Eq. (3) with respect to transverse variables $(x,y)$ is given in Ref. 19. The resultant normalized temporal correlation function of the scattered optical field is

In Eq. (7), $\theta$ is the Doppler angle [between OCT imaging direction ($z$-axis) and flow velocity vector $\mathbf{v}$]; in Eqs. (8) and (9), $\overline{A}=A+B/\rho$, $A$ and $B$ are the elements of complex valued ABCD matrix, and $A_r$, $A_i$, $B_r$, and $B_i$ are the real and imaginary parts of $A$ and $B$.

In Eq. (6), we used the Einstein–Stokes equation for spherical particle diffusivity $D_d=k_{B}T/(6\pi \eta a)$, where $k_B$ is the Boltzman constant, $T$ is the absolute temperature, $\eta$ is the liquid viscosity, and $a$ is the Brownian particle (scatterer) radius.

To interpret Eq. (5), we note that the first exponential factor is caused by the Doppler shift due to translational flow, the second exponential term is the contribution of Brownian motion, and the third term is due to the dynamic speckle fluctuations caused by translational flow. This term shows negative quadratic dependence on time in the exponent, whereas the Brownian motion term shows negative linear dependence. The physical reasons for these are the following: the scattered optical field decorrelates more rapidly the faster the particles move out of scattering volume for translational term, and move more rapidly on the wavelength scale for Brownian motion term. Particle movement is characterized by mean square displacement, which is equal to $2D\tau$ for Brownian motion along $z$-axis and to ${(v\tau )}^2$ for translational motion. Therefore, one can expect linear in time negative exponential term $\sim \exp (-2D\tau /{\lambda }^2)$ for Brownian motion and $\sim \exp [-{(v\tau /l_{sv})}^2]$ ($l_{sv}$ is characteristic dimension of scattering volume) for the translational motion term.

The flow term in Eq. (5), as seen via Eqs. (7)–(9), suggests that the correlation time due to speckle fluctuations depends on the optical system parameters (presence of matrix elements $A$ and $B$). We now evaluate the correlation function contained in Eq. (5) for a number of practical and relevant OCT sample arm geometries.

2.1.

Correlation Functions for Several Optical Coherence Tomography Sample Arm Geometries

It follows from Eqs. (5), (7), (8), and (9) that the impact of the optical system on correlation function of scattered radiation is manifest through the equivalent beam radius $w_e$. We thus now analyze several typical imaging geometries and evaluate the resulting values of $w_e$.

Table 1 shows the ABCD matrices of some optical elements²¹ typically used in the OCT sample arms. Let us now consider the particular geometries of sample arm optical system shown in Fig. 2.

Table 1

ABCD matrix of the optical elements.

Optical element	ABCD matrix
Free space: propagation over distance $z$	$M_{fs}(z)=\left(\begin{array}{cc}1& z\\ 0& 1\end{array}\right)$
Thin lens, focal length $f$	$M_{\text{lens}}(F)=\left(\begin{array}{cc}1& 0\\ -\frac{1}{F}& 1\end{array}\right)$
Gaussian diaphragm, radius $q$	$M_d(q)=\left(\begin{array}{cc}1& 0\\ -\frac{2i}{k{q}^2}& 1\end{array}\right)$
Refraction at a flat interface between media of refractive indices $n_1$ and $n_2$ (direction 1 → 2)	$M_{\mathrm{int}}(n_1,n_2)=\left(\begin{array}{cc}1& 0\\ 0& \frac{n_1}{n_2}\end{array}\right)$

Fig. 2

Geometry of the optical systems: (a) free space; (b) 4f system; (c) two-lens optical system; (d) one-lens optical system (fiber ball and/or GRIN lens); and (e) multiple lens optical system.

2.1.1.

Free space

In the case of free space shown in Fig. 2(a) [free-space collimated illumination, detection through a fiber whose end face is a distance ($L-d$) away from tissue surface], the following equation for ABCD matrix applies:

Eq. (10)

$$M=M_{fs}(L-d)M_{\mathrm{int}}(n,1)M_{fs}(dn)M_{fs}(z_1)=\left[\begin{array}{cc}1& nL-d(n-1)+z_1\\ 0& n\end{array}\right],$$

where ($L-d$) is the distance between fiber-end face (observation plane) and air–medium interface, $n$ is the medium refraction index, $d$ is the distance between position of air–medium interface and Gaussian beam waist in the air, and $z_1$ is the axial displacement of scattering volume from the beam waist.

The meaning of each matrix in Eq. (10) is explained in Table 1.

The resultant matrix elements from Eq. (10) are thus $A=1$, $B=L_1=nL-d(n-1)+z_1$, $D=n$, and $\overline{A}=1+L_1/\rho$. Plugging these into Eq. (8) yields the following equation for equivalent beam radius:

Eq. (11)

$$\frac{1}{w_e^2}=\frac{1}{w^2}+{\left(\frac{kw}{2L_1}\right)}^2{(1+\frac{L_1}{\rho })}^2.$$

The individual quantities $w$, $\rho$, and $L_1$ in Eq. (11) are $z$-dependent; however, their combination given in Eq. (11) exhibits a very weak dependence on $z$. It is not difficult to see that at the Gaussian beam waist ($z_1=0$) $w_e\approx w$. Indeed, at the beam waist the wavefront curvature $\rho =\infty$, so that the term containing it vanishes; further, in the typical OCT conditions (${\lambda }_0=1.31\mu \mathrm{m}$, $w_0=10\mu \mathrm{m}$, and $L_1\approx L=100\mathrm{mm}$), ${(k{w}_0/2L_1)}^2=0.06{\mathrm{mm}}^{-2}$, so it can be neglected when compared with $1/w_0^2={10}^4{\mathrm{mm}}^{-2}$.

2.1.2.

4f system

Figure 2(b) shows the so-called 4f system comprised two lenses. The diaphragm radius $q$ is equal to the collimated Gaussian beam radius in the plane of collimating lens, provided that Rayleigh range of the beam emerging from the fiber is much smaller than collimating lens focal distance. The beam radius in the plane of collimating lens in its turn is defined by the numerical aperture/acceptance angle of the single-mode fiber. It can also be expressed via the radius of fiber-field mode $w_f$ as ${\mathrm{NA}}_f=\lambda /(\pi w_f)$. Therefore, we obtain the following equation for the radius of (Gaussian) diaphragm:

Eq. (12)

$$q\approx {\mathrm{NA}}_f·F_c\approx \lambda F_c/(\pi w_f)=2F_c/(k{w}_f)=2F_s/(k{w}_0).$$

Here $F_c$ and $F_s$ are the focal distances of collimating and sample lenses, respectively, $w_f$ is the fiber-mode field radius, and $w_0$ is the Gaussian beam waist radius.

The ABCD matrix for backscattered radiation corresponding to this optical system can be written as

Eq. (13)

$$M=M_{fs}(F_c)M_{\text{lens}}(F_c)M_{fs}(F_c)M_d(q)M_{fs}(F_s)M_{\text{lens}}(F_s)M_{fs}(F_s-d)M_{\mathrm{int}}(n,1)M_{fs}(nd)M_{fs}(z_1)=\left(\begin{array}{cc}-\frac{F_c}{F_s}& -\frac{F_c{z}_1}{F_s}-\frac{2in{F}_c{F}_s}{k{q}^2}\\ 0& -\frac{n{F}_s}{F_c}\end{array}\right),$$

where ($F_s-d$) is the distance between sample lens and air–medium interface, $nd$ is the distance between air–medium interface and position of the Gaussian beam waist in the medium [the fact that the position of Gaussian beam waist is shifted by $(n-1)d$ in the medium relative to its position in the air follows immediately from the analysis of complex beam parameter¹⁶], and $z_1$ is the out-of-focal plane displacement of scattering volume (along the $z$-depth direction). Note that the resulting matrix does not depend on distance $d$.

By using the $A$ and $B$ matrix elements given by Eq. (13) in Eq. (8) for the equivalent beam radius (within a computational error of $\sim {10}^{-15}$), there results $w_e=w_0/\sqrt{2}$. This implies that the impact of three factors [dependence on $z$ of $w(z)$, $\rho (z)$, and matrix elements $B(z)$, $\overline{A}(z)=A+B(z)/\rho (z)$] cancels each other, producing essentially no dependence on out-of-focal plane displacement $z_1$ for equivalent beam radius $w_e$.

Although the 4f system is seldom used in the practical OCT systems, its analysis allows getting a simple solution for equivalent beam radius and, accordingly, for the temporal correlation function of scattered radiation. As shown below, this solution is actually very close to the frequently used two lens imaging system.

The analysis of other optical systems shown in Fig. 2 is performed in a similar manner. Although no simple analytical expressions are available for elements of ABCD matrices in the configurations described in Figs. 2(c)–2(e), the numerical analysis of the equivalent beam radius based on Eq. (8) is possible for these cases as well.

3.

Results and Discussion

We now compare the predictions of the developed formalism with published theoretical models.^9,10 Since OCT can yield depth-resolved measurements of correlation function and spectrum of backscattered radiation, it is useful to have model predictions not only in the focal plane of the optical system, but also at other depths as well ($\pm$ several Rayleigh ranges away from the focal plane).

Considering the results of Ref. 10 first, 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_t$. The transverse part of their point spread function is of the form $\exp [-2h_t^2(x^2+y^2)]=\exp [-2(x^2+y^2)/w^2]$. Thus, the definitions of $w$ in Ref. 10 and that in our paper differ by $\sqrt{2}$, since the transverse part of PSF in our paper has the form $\exp [-(x^2+y^2)/w^2]$ as per Eq. (1). Therefore, although the results for correlation time of scattered radiation in Ref. 10 look identical to our model for a 4f system at the axial position $z_1=0$ (i.e., in the focal plane of optical system at the beam waist), they differ by a factor of $\sqrt{2}$. For out-of-plane predictions, the difference between our model for free space and that of Ref. 10 is caused by the impact of wavefront curvature $\rho (z)$, which is not taken into account in Ref. 10. At $z=0$, $\rho (z)=\infty$, the term containing $\rho (z)$ in Eq. (1) vanishes and the two models indeed converge at the focal plane. This is illustrated in Fig. 3 (curve 1). Using our notations for results of Ref. 10, one gets $w_e\approx w_0$ with $z_1$ outside of focal plane; the model in Ref. 10 predicts $z$-dependence of $w_e$ corresponding to Gaussian beam radius $w(z)$. As noted above, our model for free space predicts essentially no $z$-dependence of the equivalent beam radius (curve 3).

Fig. 3

Comparison of published models (curves 1 and 2) for equivalent beam radius with the predictions of the model developed in this paper (curves 3–7; see Fig. 2 for relevant optical geometries). 1—the model of Ref. 10; 2—the model of Ref. 9; 3—free-space prediction of Eq. (11) of this paper; 4—one-lens imaging geometry for lensed fiber; 5—multiple lens optical system; 6–4f geometry; and 7—two-lens optical system.

Consider next the model developed in Ref. 9. It is also based on free-space solution of scattering problem, where wavefront curvature of incident beam is ignored as well. The integral equation [Eq. (12)] in Ref. 9, when evaluated in the absence of a flow velocity gradient, yields a simple analytical equation that closely resembles Eqs. (5)–(7) of this paper, with $w_e=w/\sqrt{2}$. The factor $\sqrt{2}$ in Eq. (12) of Ref. 9 is introduced empirically to account for fiber coupling efficiency. The solution thus obtained at the focal plane of optical system corresponds to our result for 4f system. Similar to Ref. 10, this model predicts Gaussian beam radius dependence of $w_e(z)=w(z)/\sqrt{2}$; it is shown as curve 2 in Fig. 3. As before, curve 3 shows virtual $z$-independence of equivalent beam radius in our model of free space given in Eq. (11). To complete the model predictions for the selected OCT sample arm geometries in Fig. 2, curve 4 shows $z$-dependence of $w_e(z)$ for the one-lens optical system formed by the lensed fiber, curve 5 corresponds to the multiple lens optical system, curve 6 corresponds to the 4f geometry, and curve 7 is for the two-lens optical system. The difference in equivalent beam radius $w_e(z)$ between these latter three optical geometries (curves 5–7) is less than 1% and therefore these curves look indistinguishable in Fig. 3. Equation (8) was used in the calculations, with parameters of optical systems given in Table 2. Thus, we predict essentially $z$-independence as one moves above or below the OCT focal plane (possible exception is the slight dependence exhibited by the lensed fiber system of curve 4). This is in sharp contrast with predictions of published models^9,10 that suggest a significant $z$-dependence shown by curves 1 and 2. Further, the actual beam radius values are off by the factor of $\sqrt{2}$, something that is accounted for naturally in our model but is introduced empirically in Ref. 9.

Table 2

Parameters of optical systems used in calculations in Fig. 3.

Let us now compare the experimental data published in Ref. 22 to the predictions of our theoretical model. It is well known that the signal component of OCT photodetector current is proportional to $2\mathrm{Re}(E_s{E}_r^{*})$, where $E_s$ and $E_r$ are the backscattered optical and reference beam fields, respectively. Since $E_r=\mathrm{const}$, the correlation function of photocurrent is directly proportional to the correlation function of its AC component.

In Ref. 22, the power spectrum of scattered radiation was studied experimentally as a function of out-of-focal plane displacement. The Gaussian beam radius $w$ was used as the fitting parameter in the theoretical model developed in Ref. 9. Figure 4 shows the values of equivalent beam radius $w_{ei}=w_i/\sqrt{2}$ as a function of out-of-focal plane (beam waist) displacement; $w_i$ are the experimentally obtained values of Gaussian beam radius as a fitting parameter.²² The dashed line shows the value of equivalent beam radius predicted by our model: $w_e=w_0/\sqrt{2}$, where $w_0=11.8\mu \mathrm{m}$ is the independently measured Gaussian beam waist radius. As seen, the predicted theoretical dependence is in good agreement with the experimental data, lending some credence to the developed formalism. The slight discrepancy (theory somewhat lower than measured data) can be explained by speckle averaging effect over the receiving aperture of the fiber-end face. Conversely, the $w(z)$ dependence of the models in Refs. 9 and 10 does not agree with the experimental data.

Fig. 4

Equivalent beam radius as a function of out-of-focal plane displacement. Circles—experimental data $w_{ei}=w_i/\sqrt{2}$ from Ref. 22; and dashed line—$w_e=w_0/\sqrt{2}$.

4.

Conclusions

We have developed a mathematical model based on the use of ABCD matrices for correlation function of optical field scattered by flowing Brownian particles in OCT conditions. The model takes into account the parameters of optical system used to collect scattered radiation onto the receiving aperture. The impact of out-of-focal plane displacement of scattering volume is considered in detail for a number of optical systems typically used in OCT sample arm configurations.

It is shown that the analytical expression for correlation function, which includes parameters of the optical system, has the same structure as that in published models based on free-space geometry, where sample arm optics are not taken into account. When other OCT sample arm configurations are considered, the resulting difference can be described by the changes in the equivalent beam radius. Several typically used OCT optical systems have been considered in detail, demonstrating that the focal plane statistics stay essentially the same throughout the OCT imaging depth. This differs significantly from the pronounced $z$-dependence predicted by existing theories and agrees well with published experimental data.