Predicting the propagation of partially coherent light in complex media is critical to optimizing existing imaging and remote sensing systems, as well as developing new measurement modalities. Such systems include confocal¹ and optical coherence tomographic² (OCT) viewing of biological tissue, optical communication through the turbulent atmosphere, and underwater imaging in littoral environments. Coherence plays an important role in these systems' performance, affecting how well incident light can be focused on a localized region, how the medium between object and pupil affects propagation, and how much scattered light can be collected from the object volume.

Traditionally, the propagation of light in complex systems (such as multiple scattering media) has been treated in one of two ways: physical optics calculations or Monte Carlo solutions to the radiative transport equation. While straightforward, the first approach is really only useful for simple geometries and weakly scattering media. On the other hand, Monte Carlo methods, while robust and able to handle more complex geometries, make simplifying assumptions and only describe incoherent propagation.

Our approach combines the physical optics description of light with stochastic methods. The resulting stochastic Huygens technique preserves the wave model of light that is essential for describing coherence phenomena, yet retains the simplicity and flexibility of Monte Carlo methods. It allows us to study complex phenomena that cannot be treated by either approach individually. This hybrid approach will facilitate the transition from empirically based design of remote sensing devices and algorithms to simulation-intensive, scientifically based concept development. We expect this formalism to enhance understanding of how light interacts with media, and to lead to better instrumentation and observational procedures.

We initially focused on problems for which we could validate our stochastic Huygens method with physical optics calculations. One such problem is the focusing of partially coherent light by an f/1 cylindrical optical system (see Figure 1). We described coherence in terms of a series of random realizations of a Gaussian-Schell line source with uniform amplitude and correlation length σ_g.

Figure 1. Illustration of f/1 focusing system.

We generated the source fields using the concept of a copula.³ Next, we employed a conventional fast Fourier transform (FFT) technique with uniformly distributed phases to generate a low-pass filtered field-realization cube.³ We then produced evolution between two statistically independent limits using a Gaussian copula that allows us to generate realizations of a random variable with arbitrary prescribed correlation (see Figure 2). The field realization cube has two spatial dimensions and a single temporal one. One can re-interpret the two spatial dimensions as indices on the sample function realizations and the temporal dimension as a spatial one instead (see video).⁴ Statistics (intensity and correlation) of this ensemble can also be viewed as a short video.⁵

Figure 2. Typical field realization cube (only real part shown here).

Figure 3 shows the concept of stochastic Huygens propagation.⁶ From randomly selected points in the source plane, rays are launched towards a second one that may contain a scattering screen, detector, or limiting aperture. At this plane, a new family of daughter rays is launched towards the next one, which may contain another limiting aperture, etc. In the final detector plane, the rays are collected at a series of discrete pixels. Note that this spatial discretization takes place only at the detector plane. Thus, this sidesteps the issue of adequate spatial sampling in any intermediate plane that occurs in physical optics calculations.

Figure 3. Illustration of stochastic sampling of Huygens wavelets.

Each of the random realizations can be propagated from source to detector. Alternatively, we can use the stochastic Huygens method to estimate the Green's function⁷ of the composite system. We can then propagate each field realization by a simple matrix multiplication. Further, if only second-order statistical descriptions of the detector field are needed, then we can propagate the cross-spectral density function through simple right and left multiplications by the Green's function. Video animations^8,9 show the field statistics (intensity, S, and coherence, μ) for the focal plane (transverse behavior) and along the axis in the vicinity of the focal plane, for varying degrees of source coherence. There is excellent agreement between the stochastic Huygens and physical optics calculations. An interesting feature of these results is that, depending on the correlation properties of the source, the coherence region at the focus can be larger or smaller than the intensity envelope.

To date, we have demonstrated that our method is fully capable of describing the effects of diffraction and partial spatial coherence. This approach has the advantage that it is able to mimic complex propagation phenomena, including difficult boundary conditions, using simple computational rules. Because the computations are based on simple rules, but are highly repetitive, they lend themselves to parallelization. We are, therefore, exploring the use of dedicated graphics processing units for accelerating the calculations. In the future, we plan to generalize our algorithms to full 3D calculations to describe polarization effects. We also plan to develop structured medium models that incorporate specified first-order properties (ray deflection, birefringence) along with second-order ones reflecting the medium's structural organization.¹⁰

We believe that this formalism for treating the propagation of light, including its coherence properties, will be useful for a diverse set of applications, particularly imaging biological tissues. Moreover, it will lead to a deeper understanding of the nature of light-matter interactions in difficult modalities such as confocal microscopy and OCT.

This work was sponsored in part by National Institutes of Health (NIH) grants CA103824 and 1R21DE016758-01A2.