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.7 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.8 Indeed, it was demonstrated recently that the determination of flow velocity and diffusivity of Brownian particles may be possible with OCT.9
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,11 and its statistical properties are different in the free space and at various planes of an imaging optical system.12
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 papers13,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 need15 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 lenses12).
An ABCD matrix approach has been able to derive solution for spatiotemporal speckle statistics for any optical system in the case of surface scattering.16 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.
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:16 The last exponential term in Eq. (1) is due to OCT coherent gating effect, where is the coherence length of the illuminating laser and is the position of scattering volume relative to beam waist. The more familiar Gaussian form of this term is an approximation. In fact, it can be viewed as a zeroth-order term of Hermite functions expansion17 of Fourier transform of source laser spectral density. The modulus of Eq. (1) gives the shape of OCT scattering volume or point spread function.
As shown in Ref. 15, if the size of scattering volume in each direction is considerably larger than 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:18
The contribution of flow can be obtained by generalizing the known expression for correlation function for surface scattering to the volumetric scattering case:19
In the general case of an ABCD optical system (paraxial approximation), the Green’s function is20
The integration of Eq. (3) with respect to axial variable can be performed easily if the variation scale of matrix elements , , and in is much larger than the source coherence length . The estimate of this scale is given in Appendix. From there, calculations show that under typical OCT conditions, the scale of matrix element -variation is , where is the Rayleigh range of laser beam in the medium. With typical values of medium refractive index , wavelength , and Gaussian beam radius at the waist , we get , whereas typical OCT coherence lengths are 5 to . The numerical estimates further show that the matrix elements and have a weak dependence on , meaning that they can be treated as constants during the integration of Eq. (3) with respect to as well.
The optical field in Eq. (1) contains other -dependent terms, specifically , , and . As is known,16 the scale of variation of these terms is ; thus all -dependent terms in Eq. (3) [apart from coherence gate and terms in the incident optical field ] 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 is given in Ref. 19. The resultant normalized temporal correlation function of the scattered optical field is
In Eq. (7), is the Doppler angle [between OCT imaging direction (-axis) and flow velocity vector ]; in Eqs. (8) and (9), , and are the elements of complex valued ABCD matrix, and , , , and are the real and imaginary parts of and .
In Eq. (6), we used the Einstein–Stokes equation for spherical particle diffusivity , where is the Boltzman constant, is the absolute temperature, is the liquid viscosity, and 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 for Brownian motion along -axis and to for translational motion. Therefore, one can expect linear in time negative exponential term for Brownian motion and ( 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 and ). We now evaluate the correlation function contained in Eq. (5) for a number of practical and relevant OCT sample arm geometries.
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 . We thus now analyze several typical imaging geometries and evaluate the resulting values of .
ABCD matrix of the optical elements.
|Optical element||ABCD matrix|
|Free space: propagation over distance|
|Thin lens, focal length|
|Gaussian diaphragm, radius|
|Refraction at a flat interface between media of refractive indices and (direction 1 → 2)|
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 () away from tissue surface], the following equation for ABCD matrix applies:
The meaning of each matrix in Eq. (10) is explained in Table 1.
The resultant matrix elements from Eq. (10) are thus , , , and . Plugging these into Eq. (8) yields the following equation for equivalent beam radius:
The individual quantities , , and in Eq. (11) are -dependent; however, their combination given in Eq. (11) exhibits a very weak dependence on . It is not difficult to see that at the Gaussian beam waist () . Indeed, at the beam waist the wavefront curvature , so that the term containing it vanishes; further, in the typical OCT conditions (, , and ), , so it can be neglected when compared with .
Figure 2(b) shows the so-called 4f system comprised two lenses. The diaphragm radius 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 as . Therefore, we obtain the following equation for the radius of (Gaussian) diaphragm:
The ABCD matrix for backscattered radiation corresponding to this optical system can be written as16], and is the out-of-focal plane displacement of scattering volume (along the -depth direction). Note that the resulting matrix does not depend on distance .
By using the and matrix elements given by Eq. (13) in Eq. (8) for the equivalent beam radius (within a computational error of ), there results . This implies that the impact of three factors [dependence on of , , and matrix elements , ] cancels each other, producing essentially no dependence on out-of-focal plane displacement for equivalent beam radius .
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.
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 ( 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 width of OCT transverse resolution” denoted as . The transverse part of their point spread function is of the form . Thus, the definitions of in Ref. 10 and that in our paper differ by , since the transverse part of PSF in our paper has the form 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 (i.e., in the focal plane of optical system at the beam waist), they differ by a factor of . 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 , which is not taken into account in Ref. 10. At , , the term containing 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 with outside of focal plane; the model in Ref. 10 predicts -dependence of corresponding to Gaussian beam radius . As noted above, our model for free space predicts essentially no -dependence of the equivalent beam radius (curve 3).
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 . The factor 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 ; it is shown as curve 2 in Fig. 3. As before, curve 3 shows virtual -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 -dependence of 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 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 -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 models9,10 that suggest a significant -dependence shown by curves 1 and 2. Further, the actual beam radius values are off by the factor of , something that is accounted for naturally in our model but is introduced empirically in Ref. 9.
Parameters of optical systems used in calculations in Fig. 3.
|Fiber-mode field radius wf (μm)||Lens focal distance, ball lens radius of curvature (mm)||Relevant distances (mm)||Medium refractive index n, ball lens refractive index nl|
|Free space [Fig. 2(a)]||4.6|
|4f system [Fig. 2(b)]||4.6|
|Two-lens optical system [Fig. 2(c)]||4.6|
|One-lens optical system [Fig. 2(d)]||4.6|
|Multiple lens optical system [Fig. 2(e)]||4.0|
|For all calculations, central wavelength ;|
|Gaussian beam radius at the waist (both in air and in medium)|
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 , where and are the backscattered optical and reference beam fields, respectively. Since , 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 was used as the fitting parameter in the theoretical model developed in Ref. 9. Figure 4 shows the values of equivalent beam radius as a function of out-of-focal plane (beam waist) displacement; are the experimentally obtained values of Gaussian beam radius as a fitting parameter.22 The dashed line shows the value of equivalent beam radius predicted by our model: , where 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 dependence of the models in Refs. 9 and 10 does not agree with the experimental data.
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 -dependence predicted by existing theories and agrees well with published experimental data.
Examination of the -dependence of the Quantities in Eq. (4) to Enable the Integration of Eq. (3)
Considering the 4f system first, one gets the following equations for its matrix elements from Eq. (13):
Equation (14) shows that only element depends on . Let us estimate the scale of variation for [since enters Eq. (4) for Green’s function as ].
By putting from Eq. (12) into equation for in Eq. (14), one gets
Here is the Rayleigh range of the Gaussian beam in the medium. Equation (16) shows that the -spatial scale of variation is .
The numerical estimates show that for two-lens and multiple lens optical systems, the dependence of elements and on variable is weak on the scale of coherence length and can be safely neglected. has approximately the same scale of variation as for 4f system.
A similar calculation for one-lens optical system with a ball lens spliced to the fiber shows that the spatial scale of element variation is given by
By putting Eqs. (1) and (4) into Eq. (3), we get the following equation for the axial part of the integrand in Eq. (3):
It follows from Eq. (18) that its point of extremum is at . This means that the slowly changing terms in Eq. (4) should be evaluated at this point; they are becoming time dependent.
Let us estimate the magnitude of time-dependent term . Since we are not interested in times exceeding the correlation time of OCT signal, we will estimate the value of , where is the time corresponding to level of temporal correlation function. We get from Eqs. (5) and (7): , so that . This means that the term and can be neglected safely.
Recapping, we show that the transverse part of Green’s function in Eq. (4), given by elements , , and , has a spatial scale of variation along -axis for all optical systems considered. Therefore, evaluation of integral in Eq. (3) with respect to variable in the limit of can be performed by assuming the transverse part of Green’s function is approximately constant in .
This study was supported by the Canadian Institutes of Health Research (Grant No. 126172), the Natural Sciences and Engineering Research Council of Canada (194509), and the Ministry of Education and Science of the Russian Federation (14.B25.31.0015).
Alex Vitkin is a professor of medical biophysics and radiation oncology at the University of Toronto, a senior scientist at the University Health Network, and a clinical physicist at Princess Margaret Cancer Centre (all in Toronto, Canada). He has published ~150 papers/book chapters on biophotonics, primarily on tissue polarimetry and optical coherence tomography. He is also a fellow of OSA and SPIE, and a professor at the Nizhny Novgorod State Medical Academy (Russia).