Open Access
Add to BibSonomyAbstract
Monte Carlo techniques are the gold standard for studying light propagation in turbid media. Traditional Monte Carlo techniques are unable to include wave effects, such as diffraction; thus, these methods are unsuitable for exploring focusing geometries where a significant ballistic component remains at the focal plane. Here, a method is presented for accurately simulating photon propagation at the focal plane, in the context of a traditional Monte Carlo simulation. This is accomplished by propagating ballistic photons along trajectories predicted by Gaussian optics until they undergo an initial scattering event, after which, they are propagated through the medium by a traditional Monte Carlo technique. Solving a known problem by building upon an existing Monte Carlo implementation allows this method to be easily implemented in a wide variety of existing Monte Carlo simulations, greatly improving the accuracy of those models for studying dynamics in a focusing geometry.
© 2015 Optical Society of America
1. Introduction
Monte Carlo methods have become the gold standard for analyzing light propagation in turbid media [1, 2], and have been used to investigate an enormous range of problems including thermal dynamics [3], single- and multi-photon fluorescence [4–7], spontaneous Raman scattering [8, 9], and stimulated Raman scattering [10, 11]. Monte Carlo techniques have been instrumental tools in the understanding of light propagation in turbid media which has lead to the development of deep tissue microscopy [12–15] and random lasing [16, 17].
Monte Carlo simulations approximate elastic scattering as a stochastic process in which photon packets propagate through the turbid media in a random walk. By assigning probability distribution functions to the relevant physical processes, such as absorption and scattering, Monte Carlo simulations are able to calculate many physically relevant quantities, like the intensity distribution of light in the reflection or transmission geometry. Mathematically, Monte Carlo simulations are equivalent to solving the radiative transport equation, and while Monte Carlo simulations are typically more computationally demanding than numerically solving the transport equations they are often preferred because they are capable of taking into account complicated sample geometries and boundary conditions.
One well known problem associated with traditional Monte Carlo simulations is their inability to include diffraction due to treating light as a particle [18]. This severely limits their ability to study the dynamics of focusing beams in situations where some ballistic, unscattered, photons remain at the focal plane, because these photons will focus to a point instead of a diffraction limited spot size. Nevertheless, several models have been introduced to address this problem, but such methods generally rely on adopting an initial distribution of photons that will result in a Gaussian profile at the focus [6,19–23]. Such a method can give the correct profile at the focus in the absence of scattering; however, it will, in general, give the wrong distribution everywhere else. To date, the only solutions are to use electric field Monte Carlo (EMC) [18,24] simulations or wave Monte Carlo (WMC) [25, 26] simulations which are capable of including wave effects such as diffraction; however, adoption of these methods would require significant reworking of existing traditional particle-based Monte Carlo models. Thus, a way to adapt a traditional Monte Carlo simulation to be able to accurately simulate a focusing beam is still needed.
In this paper, a simple generalization of a traditional Monte Carlo method, in which ballistic photons follow paths derived from Gaussian optics is presented. Thus, diffraction-limited spots can be accurately described within the framework of traditional Monte Carlo simulations. We have based our implementation on a traditional Monte Carlo multi-layer (MCML) [2]; however, the approach taken here is general enough that application to other Monte Carlo techniques should be straightforward.
2. Derivation
To include Gaussian optics into a standard Monte Carlo simulation we will need to modify only the initial distribution of the pulse and the behavior of the unscattered, or ballistic, photons. The Monte Carlo method presented in this work is very similar to the standard MCML model. For this reason we will describe it only briefly; for additional details please see Ref. [2]. The photons are initialized outside the sample and distributed into Gaussian spatial and temporal distributions (Note that this will be modified later to account for a focusing lens in the system). For the photons outside the sample a temporal delay, Δτ, is used such that the center of the Gaussian pulse will arrive at the surface of the sample after a time Δτ. The photons are then propagated through the sample with scattering lengths sampled from an exponential distribution and scattering angles sampled from the Henyey-Greenstein distribution. Absorption is included by assigning each photon an initial weight of unity which is depleted each step according to W = W exp(−μa * Δr), where μa is the absorption coefficient, and Δr is the distance traveled by the photon in that step. To reduce the simulation noise, photon propagation is terminated according to a Russian Roulette-type process once the photon weight falls below a certain threshold value [2].
The scalar electric field of a monochromatic, uniformly polarized Gaussian beam is given by
where k is the wavenumber, w(z) is the beam radius, R(z) is the radius of curvature, and ζ(z) is the Gouy phase. The radius of curvature of a Gaussian beam is given by Here, zf is the z location of the focus with z = 0 corresponding to the surface of the sample. We take z to be positive inside the sample and negative to the reflection side of the sample. The Rayleigh length, zR, is defined as where λ is the wavelength of light. The beam radius is given by where w0 is the beam radius at the focus. In the paraxial approximation w0 = (2/π)(λf/D), where f is the focal length of the lens used, and D is the diameter of the beam at the lens.To define the initial beam we must define both the initial velocity and the initial position of the photons. We arrive at the correct initial photon distribution by starting with a Gaussian spatio-temporal pulse in collimated space and then applying a variable transformation which transforms the pulse into focusing space by accounting for the effects of a lens. We should note that for this approach to be valid, it is assumed that the pulse is initialized far from the focus. In essence, this approach takes rays which would ordinarily run parallel to the z-axis in collimated space and makes them cross at the focus in focusing space. This process is schematically depicted in Fig. 1 and can be described as a simple coordinate change in the initial position distribution of the photons as
Here, f is the focal length of the lens and x′, y′, and z′ are the photon coordinates in collimated space, which for our purposes are generated by sampling a Gaussian distribution. By using these transformations we could, in principle, generalize the method outlined here to non-Gaussian beams in either space or time. However, generalization to non-Gaussian spatial beams would require using a different propagation method if their trajectories could not be satisfactorily described in the context of Gaussian optics.
Fig. 1 (A) Conceptual diagram of the transformation from collimated space to focusing space. (B) Schematic illustrating the interplay between the wavefront and the direction vector of the photon. The dashed lines represent the trajectories that would occur in a traditional Monte Carlo.
A ray, or photon packet in the Monte Carlo simulation, has a direction perpendicular to the radius of curvature of the beam (see Fig. 1). This direction unit vector is given by
To simulate a focusing beam we use Eq. (5) to assign the direction vector of each photon at each timestep until it undergoes its first scattering event. Once the first scattering event occurs, the photon is propagated through the remainder of the Monte Carlo simulation in the traditional manner. In our implementation we adhere to the commonly used Monte Carlo multilayer method [2]. This technique of propagating photons by Gaussian dynamics prior to the first scattering event solves the problem of ballistic photons focusing to a point, and allows focusing geometries to be accurately modeled in Monte Carlo simulations with much greater accuracy. To avoid numerical instabilities, we define T (z) = 1/R(z) instead of finding R(z) directly.To approximate the trajectories of the ballistic photons we must solve the equation of motion ṙ = (c/n)v̂. This first-order differential equation can be solved with many numerical methods. In practice, because the sign of the curvature does not change, low-order methods, such as Euler’s, perform poorly. The local error of a method such as Euler’s will always have the same sign, resulting in a large global error accrued by such a technique. Thus, we use a 4th-order Runge-Kutta method [27] in our implementation with good results. To improve the accuracy near the focus, and minimize the run time far away from the focus, we vary the step size, Δt, depending on the radius of curvature of the trajectories. To accomplish this, we allow the step size to change depending on the largest radius of curvature of the three components. We define the step size to be
where the factor of c/n has been inserted to make ε a unitless parameter. ẍ, ÿ, and z̈ can be found by straightforward differentiation of the equation of motion3. Results and Discussion
To demonstrate the ability of this method to generate the proper distribution at the focus, and to confirm that the process is numerically stable, we turn scattering and absorption off in the Monte Carlo simulation. The results of simulating a 10-mm 1/e2-diameter, 1-ps full-width at half-maximum temporal width beam being focused by a 25.4-mm lens are shown in Fig. 2. Here, the beam was given a wavelength of 1.0 μm and was represented by 106 photons. The simulation was repeated 10 times with the same parameters (but using a different seed for the random number generator) and averaged. The theoretical 1/e2-diameter at the focus in the paraxial approximation is 3.2340 μm, and the simulation gives a value of 3.2326 μm with a step size parameter, ε, of 10−6. The temporal pulse width given by the simulation was 1.0005 ps. Both of these results demonstrate clearly that this method is capable of accurately producing the correct distribution at the focus with no temporal aberrations.
Fig. 2 Spatial, (A), and temporal, (B), profile at the focus of a 1-ps, 1-μm beam focused with a 25.4-mm focal length lens with scattering turned off in the Monte Carlo simulation. The incident beam was 10.0-mm 1/e2-diameter at the lens. The green lines are Gaussian fits to the data.
To demonstrate this method’s performance in scattering media we keep the same parameters as before; however, we now propagate the beam through a 1-mm-thick turbid sample. The sample is set so that the focus of the lens corresponds to the back surface of the sample. We vary the scattering coefficient of the sample to demonstrate how the beam profile at the focus changes with scattering in the regime where the mean free path is on the order of the sample thickness. The results are shown in Fig. 3. In these simulations the scattering is assumed to be largely anisotropic with g = 〈cos(θ)〉 = 0.9. It can be readily seen that both the spatial and temporal profiles undergo significant changes in the region where the scattering coefficient is on the order of the sample thickness. When only a traditional Monte Carlo approach is used without keeping track of the Gaussian trajectories, the spatial profiles predicted are completely wrong, but the temporal profiles are still accurately reproduced because the overall pathlength between Gaussian trajectories and ray optics are very similar.
Fig. 3 Spatial, (A) and (C), and temporal, (B) and (D), profiles at the focus of 1 ps, 1 μm beams focused with a 25.4 mm focal length lens through a 1 mm thick index matched sample with the scattering coefficient given. In each case, g = 0.9 and the medium is assumed to have zero absorption. The focal plane coincides with the back plane of the medium in the Monte Carlo simulation. (A) and (B) track the Gaussian trajectories while (C) and (D) use a traditional Monte Carlo approach. Note that the tiny spec in the middle of each panel of (C) is the distribution. The spatial profiles in (A) have been normalized so that the average value of each plot is equal. (C) is scaled by the same factor as those in (A) to illustrate that traditional Monte Carlo dramatically over predicts the focal plane intensity. The scale of the inset in (C) goes from −0.2 μm to 0.2 μm.
The examples given in this paper are focused on the regime where the reduced scattering coefficient, μ′s = (1 − g)μs, is smaller but comparable to the inverse of the focal depth, L. This regime is where traditional Monte Carlo approaches fail to give accurate results. In this regime, scattering is not large enough to sufficiently diffuse the pulse before the focus, thus propagating the pulse using geometrical optics will result in the residual ballistic light focusing to a point, causing the simulation to greatly overestimate the intensity of the pulse at the focus. As scattering is increased beyond this regime, the focus will be sufficiently degraded by multiple scattering and the errors from geometrical optics will begin to diminish, but this will not typically happen until the pulse has sufficiently entered the diffusive regime (μ′sL ≫ 1).
All of the results outlined here assume an index-matched sample; however, this method could be generalized to include the effects of different index layers by employing Snell’s law at the interface for the rays [28]. The effects of refraction would lead to new parameters, zR and zf, that would describe the new beam. If one would like to include aberrations that would result from the interface each photon could be assigned a different zR and zf.
Due to the stiffness of the equations for the trajectories, small time-steps generally must be used if accurate results for the spatial profiles are required. This can lead to the technique being relatively computationally intensive, especially where the curvature is largest near the focal plane. In our implementation, like many modern Monte Carlo implementations, we make use of the massively parallel architecture of modern graphics processing units (GPUs) [29, 30]. By utilizing the power of modern GPU’s these simulations generally take anywhere from a few seconds to a few minutes to simulate a single pulse consisting of 106 photons, depending on the parameters used.
4. Conclusion
A technique for accurately simulating a focusing beam in a traditional Monte Carlo simulation has been presented. This technique is simple to incorporate within an existing Monte Carlo simulation, and is capable of accurately simulating focusing Gaussian beams in a turbid medium with any amount of scattering. It will be immediately useful in a wide variety of applications ranging from understanding the safety of focusing pulses in tissue to optimizing the depth and resolution of deep-tissue microscopy techniques. The source code for this paper will be posted to https://github.com/hokrbh/MCFocusing.git once it is cleared for public release.
Acknowledgments
This work was partially supported by National Science Foundation Grants ECCS-1250360, DBI-1250361, CBET-1250363, PHY-1241032 (INSPIRE CREATIV), PHY-1068554, EEC-0540832 (MIRTHE ERC) and the Robert A. Welch Foundation ( Award A-1261). BHH would like to acknowledge a graduate fellowship from the Department of Defense Science, Mathematics and Research for Transformation (SMART) fellowship program. JNB and RJT would like to acknowledge support from AFRL 711 HPW. Research performed by TASC Inc. and Nanohmics Inc. was conducted under USAF Contract Number FA8650-14-D-6519.
References and links
1. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10, 824–830 (1983). [CrossRef] [PubMed]
2. L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]
3. J. N. Bixler, B. H. Hokr, M. L. Denton, G. D. Noojin, A. D. Shingledecker, H. T. Beier, R. J. Thomas, B. A. Rockwell, and V. V. Yakovlev, “Assessment of tissue heating under tunable near-infrared radiation.” J. Biomed. Opt. 19, 070501 (2014). [CrossRef]
4. J. Y. Qu, C. E. MacAulay, S. Lam, and B. Palcic, “Laser-induced fluorescence spectroscopy at endoscopy: tissue optics, Monte Carlo modeling, and in vivo measurements,” Opt. Eng. 34, 3334–3343 (1995). [CrossRef]
5. Q. Liu, C. Zhu, and N. Ramanujam, “Experimental validation of Monte Carlo modeling of fluorescence in tissues in the UV-visible spectrum,” J. Biomed. Opt. 8, 223–236 (2003). [CrossRef] [PubMed]
6. C. M. Blanca and C. Saloma, “Monte Carlo analysis of two-photon fluorescence imaging through a scattering medium,” Appl. Opt. 37, 8092–8102 (1998). [CrossRef]
7. X. Deng and M. Gu, “Penetration depth of single-, two-, and three-photon fluorescence microscopic imaging through human cortex structures: Monte Carlo simulation,” Appl. Opt. 42, 3321–3329 (2003). [CrossRef] [PubMed]
8. N. Everall, T. Hahn, P. Matousek, A. W. Parker, and M. Towrie, “Photon migration in Raman spectroscopy,” Appl. Spectrosc. 58, 591–597 (2004). [CrossRef] [PubMed]
9. B. H. Hokr and V. V. Yakovlev, “Raman signal enhancement via elastic light scattering,” Opt. Express 21, 11757–11762 (2013). [CrossRef] [PubMed]
10. B. H. Hokr, J. N. Bixler, and V. V. Yakovlev, “Higher order processes in random Raman lasing,” Appl. Phys. A 117, 681–685 (2014). [CrossRef]
11. B. H. Hokr, V. V. Yakovlev, and M. O. Scully, “Efficient time-dependent Monte Carlo simulations of stimulated raman scattering in a turbid medium,” ACS Photonics 1, 1322–1329 (2014). [CrossRef]
12. I V Meglinskii, A N Bashkatov, E A Genina, D Yu Churmakov, and V. V. Tuchin, “Study of the possibility of increasing the probing depth by the method of reflection confocal microscopy upon immersion clearing of near-surface human skin layers,” Quantum Electron. 32, 875–882 (2002). [CrossRef]
13. F. Helmchen and W. Denk, “Deep tissue two-photon microscopy,” Nat. Methods 2, 932–940 (2005). [CrossRef] [PubMed]
14. K. Welsher, S. P. Sherlock, and H. Dai, “Deep-tissue anatomical imaging of mice using carbon nanotube fluorophores in the second near-infrared window,” Proc. Natl. Acad. Sci. 108, 8943–8948 (2011). [CrossRef] [PubMed]
15. N. G. Horton, K. Wang, D. Kobat, C. G. Clark, F. W. Wise, C. B. Schaffer, and C. Xu, “In vivo three-photon microscopy of subcortical structures within an intact mouse brain,” Nat. Photonics 7, 205–209 (2013). [CrossRef]
16. H. Cao, Y. G. Zhao, S. T. Ho, E. W. Seelig, Q. H. Wang, and R. P. H. Chang, “Random laser action in semiconductor powder,” Phys. Rev. Lett. 82, 2278–2281 (1999). [CrossRef]
17. B. H. Hokr, J. N. Bixler, M. Cone, J. D. Mason, H. T. Beier, G. D. Noojin, G. I. Petrov, L. A. Golovan, R. J. Thomas, B. A. Rockwell, and V. V. Yakovlev, “Bright emission from a random Raman laser,” Nat. Commun. 5, 4356 (2014). [CrossRef] [PubMed]
18. C. K. Hayakawa, E. O. Potma, and V. Venugopalan, “Electric field Monte Carlo simulations of focal field distributions produced by tightly focused laser beams in tissues,” Biomed. Opt. Express 2, 278–299 (2011). [CrossRef] [PubMed]
19. J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A, Opt. Image Sci. 13, 952–961 (1996). [CrossRef]
20. Z. Song, K. Dong, X. H. Hu, and J. Q. Lu, “Monte Carlo Simulation of Converging Laser Beams Propagating in Biological Materials,” Appl. Opt. 38, 2944 (1999). [CrossRef]
21. X. Gan and M. Gu, “Effective point-spread function for fast image modeling and processing in microscopic imaging through turbid media,” Opt. Lett. 24, 741 (1999). [CrossRef]
22. X. Deng, X. Wang, H. Liu, Z. Zhuang, and Z. Guo, “Simulation study of second-harmonic microscopic imaging signals through tissue-like turbid media.” J. Biomed. Opt. 11, 024013 (2006). [CrossRef]
23. A. Leray, C. Odin, E. Huguet, F. Amblard, and Y. Le Grand, “Spatially distributed two-photon excitation fluorescence in scattering media: experiments and time-resolved Monte Carlo simulations,” Opt. Commun. 272, 269–278 (2007). [CrossRef]
24. F. Cai and S. He, “Electric field Monte Carlo simulation of focused stimulated emission depletion beam, radially and azimuthally polarized beams for in vivo deep bioimaging.” J. Biomed. Opt. 19, 11022 (2014). [CrossRef]
25. V. P. Kandidov, “Monte Carlo method in nonlinear statistical optics,” Physics-Uspekhi 39, 1243–1272 (1996). [CrossRef]
26. V P Kandidov, V O Militsin, A. V. Bykov, and A. V P., “Application of corpuscular and wave Monte-Carlo methods in optics of dispersive media,” Quantum Electron. 36, 1003–1008 (2006). [CrossRef]
27. R. L. Burden and J. D. Faires, Numerical Analysis (Brookes/Cole, Boston, MA, 2010), 9th ed.
28. P. W. Milonni and J. H. Eberly, Lasers (Wiley, 1988), 1st ed.
29. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178 (2009). [CrossRef] [PubMed]
30. A. Doronin and I. Meglinski, “Online object oriented Monte Carlo computational tool for the needs of biomedical optics,” Biomed. Opt. Express 2, 2461–2469 (2011). [CrossRef] [PubMed]
