Electric field Monte Carlo simulations of focal field distributions produced by tightly focused laser beams in tissues

Carole K. Hayakawa; Eric O. Potma; Vasan Venugopalan

doi:10.1364/BOE.2.000278

1. Introduction

The image contrast in laser scanning microscopy of biological samples depends critically on the ability to form a tightly focused spot in the specimen. The presence of scattering and absorption in biological materials generally affects the laser focus quality [1]. Scattering, in particular, is significant in biological tissues, and decreases the amount of radiation available for the formation of the focal spot at greater depths. Moreover, tissue scattering introduces both spatial distortions of the focal volume and a depolarization of the incident light. While some of the negative effects of tissue scattering can be mitigated by using higher incident light intensities [2], longer excitation wavelengths [3], or adaptive optics techniques [4, 5], the understanding of the fundamental mechanisms that link tissue properties and focus quality can offer clues towards the development of better image correction schemes.

A predictive model that connects the electric field characteristics in the focal volume to tissue parameters would be quite valuable. Generally, light propagation in complex media involves solving Maxwell’s equations of electromagnetic radiation with appropriate boundary conditions. This approach is taken when applying perturbation theory to calculate the scattered electric field [6], or when numerically solving Maxwell’s equations using finite-difference time domain (FDTD) methods [7–9]. FDTD methods have been applied recently to focused beam geometries [10], and used to study focal field distortions introduced by cellular structures [11]. However, direct solutions for the electric field is very computationally intensive when applied to tissue scattering problems, and such calculations are not easily implemented beyond specific deterministic structures of the scattering material.

An alternate approach is found in the radiative transfer equation (RTE), which models the incoherent propagation of light through scattering media [12, 13]. The RTE can be solved by direct numerical integration [14] or by Monte Carlo based methods [15–18]. The Monte Carlo (MC) method is particularly attractive, as it allows the simulation of light propagation through the tissue using probability density functions (pdf) that govern the probability that a photon interacts with the tissue as it propagates. These pdfs are parameterized by experimentally-accessible tissue parameters such as the scattering coefficient µ_s, the absorption coefficient µ_a and the single-scattering phase function p(θ). The MC method has been very successful at predicting light propagation in tissues based on general material parameters, in particular in the transport and diffusive regimes of propagation [19]. Nonetheless, because the wave character of light is generally ignored in this approach, diffraction of light is not included in MC models and the simulation of the diffraction-limited focal volume is intrinsically problematic. Several models have been introduced that provide effective solutions to this problem. Such models generally adopt initial distribution functions for the photon particles that mimic the light distribution of a Gaussian-shaped focal volume, while the light propagation still proceeds in an incoherent manner [20–29]. In addition, some MC studies have incorporated effective phase retardation functions in focused light geometries to calculate the loss of interference efficiency in optical coherence tomography (OCT) [30–33].

While MC models that incorporate effective focal volume geometries have reproduced some experimentally-observed trends, they fail to make a general connection between tissue parameters and the electric field characteristics in the vicinity of the focal volume, including its amplitude, phase, and polarization state. The characterization of the focal volume in terms of the electric field, as opposed to incoherent photon particles, is crucial for modeling the imaging properties of coherent imaging techniques such as OCT, harmonic generation microscopy, and coherent Raman microscopy. Several MC techniques have been developed that model the amplitude and phase of the electric field as it propagates through the medium. Fisher and coworkers modeled light propagation in terms of Huygens wavelets, which evolve through MC sampling [34]. Another approach is based on the decomposition of wavefronts into an angular spectrum of plane waves, which are subsequently propagated in a Monte Carlo fashion. Daria and co-workers have used such an approach to study light propagation in tissues [35]. A more formal implementation of this technique was developed by Xu, who applied the plane wave electric field Monte Carlo (EMC) method to study the spatial coherence of light in backscattered geometries [36–38].

The angular spectrum representation of diffraction theory can be integrated with the plane wave EMC method to simulate the propagation of a focused wavefront in scattering media. This is possible because the amplitude, phase, and polarization state are retained in the plane wave EMC model which allows for diffraction effects to be included in the description of the focal fields. This hybrid approach has recently been used to study the amplitude and phase of focal fields as a function of depth in scattering media [39]. In this work, we apply this method to establish general trends between experimentally accessible tissue parameters and the amplitude loss, spatial distortion, and polarization loss of focal fields. Within the random scattering approximation, we examine the effects of µ_s and the scattering anisotropy g on the quality of the focal volume as a function of focusing depth and the numerical aperture (NA) of the lens.

2. Theory

We use angular spectrum representation of diffraction theory in combination with electric field Monte Carlo (EMC) to model the formation of the focal volume in the sample. Scattering effects in the turbid medium are implemented by propagating the plane waves of the angular spectrum through an EMC simulation. The EMC simulation accounts for amplitude, phase, and polarization state changes of the electric field introduced by scattering and absorption events in the medium.

2.1. Angular spectrum representation of diffraction theory

In the angular spectrum representation the wavefront is decomposed into a spectrum of plane waves, each of which is characterized by a wave vector k. The angular spectrum representation of E in the vicinity of the focal volume is [40, 41]:

Fig. 1. Schematic of the diffraction geometry. The wavefront of the initial field E _far is modified to E^d _far, which captures the effects of a given medium. Waves launched from a Lambertian source (symbolized by semi-circle) with a wave vector k are allowed to scatter in a medium of thickness T, and the amplitude and phase at each exit wave vector k′ is determined.

Download Full Size | PDF

E_{f} (x, y, z) = \frac{i f e^{- i k f}}{2 π} \int \int_{(k_{x}^{2} + k_{y}^{2}) \leq k^{2}} E_{far}^{d} (k_{x}^{'}, k_{y}^{'}) e^{i (k_{x}^{'} x + k_{y}^{'} y + k_{z}^{'} z)} \frac{1}{k_{z}^{'}} d k_{x}^{'} d k_{y}^{'}

where f is focal length of the lens and E ^d _far is the refracted field at the lens surface. In cylindrical coordinates for the focal field, Eq. (1) is written as

E_{f} (ρ, φ, z) = \frac{i k f e^{- i k f}}{2 π} \int_{ϕ = 0}^{2 π} \int_{θ = 0}^{θ_{max}} E_{far}^{d} (θ', ϕ') e^{i k z \cos θ'} e^{i k ρ \sin θ' \cos (ϕ' - φ)} \sin θ' d θ' d ϕ' .

We incorporate turbidity by introducing a response function that represents the amplitude decay and phase delay:

E_{far}^{d} (θ', ϕ') = \int_{ϕ = 0}^{2 π} \int_{θ = 0}^{\frac{π}{2}} G (θ, ϕ \to θ', ϕ') E_{far} (θ, ϕ) \sin θ d θ d ϕ

where G(θ, ϕ →θ′, ϕ′) is called the coherent angular dispersion function (CADF), and E _far is the unperturbed field at the lens surface:

E_{far} (θ, ϕ) = {(\frac{n_{1}}{n_{2}})}^{\frac{1}{2}} \sqrt{\cos θ} (\begin{matrix} \cos ϕ \cos θ \cos (ϕ - γ) + \sin ϕ \sin (ϕ - γ) \\ \sin ϕ \cos θ \cos (ϕ - γ) - \cos ϕ \sin (ϕ - γ) \\ \sin θ \cos (ϕ - γ) \end{matrix}) ∣ E_{inc} (θ, ϕ) ∣ .

Here the unrefracted field incident at the lens aperture is written as ∣E _inc(θ, ϕ)∣.

The EMC determines G(θ, ϕ → θ′,ϕ′), that is, the amplitude loss and phase retardation associated with the scattering of an incident wave vector k to an exiting wave vector k′ (see Fig. 1). For a transparent sample, no scattering would occur and thus waves incident in k would exit at k′ = k without attenuation resulting in a G equivalent to the identity matrix.

For turbid samples, scattering produces nonzero off diagonal elements in the CADF. The resulting G matrix describes how the incident wavefront is altered due to the effects of scattering. This response function is then inserted into the diffraction equation, Eq. (3), to provide a full description of the resulting electric field that incorporates the effects of sample turbidity.

2.2. Monte Carlo simulations

We use an electric field Monte Carlo simulation to determine the CADF G(θ, ϕ → θ′,ϕ′). A Lambertian source launches plane waves at the slab surface and the initial wave vector k = (k_x, k_y) is noted. The initial local coordinate system is described by (m̂, n̂, ŝ) where ŝ is the propagation direction of the plane wave and m̂ and n̂ are unit vectors collinear with the parallel and perpendicular components of the electric field, respectively. We set the incident electric field to be linearly polarized along m̂, i.e., E = E _∥ m̂ + E _⊥ n̂ where E _∥ = 1+0i and E _⊥ = 0+0i and set the wave weight to W = 1. Each plane wave with wave vector k is launched and allowed to propagate and scatter in the medium. After each scattering/absorption event, the coordinate system is updated according to [36]:

(\begin{matrix} \hat{m}' \\ \hat{n}' \\ \hat{s}' \end{matrix}) = M (θ, ϕ) (\begin{matrix} \hat{m} \\ \hat{n} \\ \hat{s} \end{matrix})

where the coordinate transformation is given as:

M (θ, ϕ) = (\begin{matrix} \cos θ \cos ϕ & \cos θ \sin ϕ & - \sin ϕ \\ - \sin ϕ & \cos ϕ & 0 \\ \sin θ \cos ϕ & \sin θ \sin ϕ & \cos θ \end{matrix}) .

Here θ is the scattering angle and ϕ the azimuthal angle. Interactions between the plane waves and the turbid medium has the effect of altering the propagation direction ŝ of the plane wave and its corresponding projections on m̂ and n̂.

Intercollision distances Discrete absorption weighting [42] is used to model absorption which determines intercollision distances based on exponential distribution a mean length 1/µ_t and attenuates the wave at each interaction by (µ_s/µ_t), where µ_t = µ_a+µ_s.

Scattering angles The scattering angles (θ, ϕ) are determined from a joint distribution phase function p(θ, ϕ) using a method from Xu [36]:

p (θ, ϕ) = \frac{F (θ, ϕ)}{π x^{2} Q_{sca}}

where Q _sca is determined from Mie scattering calculations [43] and x is the size parameter ≡ (2πna/λ), with n as the refractive index of the medium and a as the particle radius. In Eq. (9), F(θ, ϕ) is given by

F (θ, ϕ) = ({∣ S_{2} ∣}^{2} \cos^{2} ϕ + {∣ S_{1} ∣}^{2} \sin^{2} ϕ) {∣ E_{∥} ∣}^{2} + ({∣ S_{2} ∣}^{2} \sin^{2} ϕ + {∣ S_{1} ∣}^{2} \cos^{2} ϕ) {∣ E_{⊥} ∣}^{2} +

where S ₁(θ) and S ₂(θ) are defined by Mie scattering calculations. The scattering angle θ is sampled from

p (θ) = \int_{0}^{2 π} p (θ, ϕ) d ϕ = \frac{{∣ S_{1} (θ) ∣}^{2} + {∣ S_{2} (θ) ∣}^{2}}{x^{2} Q_{sca}} .

The azimuthal angle ϕ is determined using rejection sampling from the conditional probability p(ϕ∣θ) = p(θ, ϕ)/p(θ) where p(θ) given by Eq. (9). After each scattering event, the coordinate system is updated according to Eq. (5) and the electric field is updated according to:

(\begin{matrix} E_{∥}^{'} \\ E_{⊥}^{'} \end{matrix}) = L (θ, ϕ) (\begin{matrix} E_{∥} \\ E_{⊥} \end{matrix})

with

L (θ, ϕ) = \frac{1}{\sqrt{F (θ, ϕ)}} (\begin{matrix} S_{2} (θ) \cos ϕ & S_{2} (θ) \sin ϕ \\ - S_{1} (θ) \sin ϕ & S_{1} (θ) \cos ϕ \end{matrix}) .

Tallies The propagation of the wave continues until it exits the top of the slab at the focal plane. Upon crossing the focal plane, we determine the exit angle with respect to θ subdivisions or bins in [0,π/2]. The resulting wave weight W_j for each wave j is added to the appropriate bin to determine the wave tally for bin p. The phase delay is determined from path length information. At the exit, we convert the path length into time using, t_j =(d_j−d_b)/(c/n), where d_b is the “ballistic” distance from the wavefront launch position/angle to the focal plane, d_j is the actual wavefront path length upon arrival at the focal plane, c is the speed of light, and n is the refractive index of medium. To calculate the phase delay we define t _cycle = λn/c to denote the time of a full cycle (2π). Each wave that enters angular bin p at time t_j has a phase ϕ_j

ϕ_{j} = [(t_{j} / t_{cycle}) - floor (t_{j} / t_{cycle})] * 2 π .

To determine the change to the electric field coordinate system, we must determine the change to the exiting wave coordinates relative to the initial coordinates. Let (m̂, n̂, ŝ) denote the initial electric field coordinates and (m̂ _f, n̂ _f, ŝ _f) the final wave coordinates. After n scattering events, the local coordinate system is:

(\begin{matrix} {\hat{m}}_{f} \\ {\hat{n}}_{f} \\ {\hat{s}}_{f} \end{matrix}) = M (θ_{p}, ϕ_{p}) M^{f} (\begin{matrix} \hat{m} \\ \hat{n} \\ \hat{s} \end{matrix})

with M ^f = Πⁿ _i=1 M(θ_i,ϕ_i) and M(θ_p,ϕ_p) representing the transformation back onto the original coordinate system with rotation angles θ_p = arctan(M^f ₁₃/M^f ₃₃) and ϕ_p = arctan(M^f ₂₃/M^f ₁₃). The final electric field, E ^f = E^f _∥ m̂ _f + E^f _⊥ n̂ _f, is then written as:

(\begin{matrix} E_{∥}^{f} \\ E_{⊥}^{f} \end{matrix}) = (\begin{matrix} \cos ϕ_{p} & \sin ϕ_{p} \\ - \sin ϕ_{p} & \cos ϕ_{p} \end{matrix}) L^{f} (\begin{matrix} E_{∥} \\ E_{⊥} \end{matrix})

where L ^f = Πⁿ _i=1 L(θ_i,ϕ_i) contains all the scattering-induced coordinate transformations. Note that in the EMC simulations, the electric field components E _∥ and E _⊥ are complex and are thus characterized by both an amplitude and a phase. For each wave, the phase change is determined at each scattering event and the propagation phase in between scattering events is calculated, from which the final phase of the wave upon exiting the slab is determined. In this fashion, the amplitude and phase wavefront can be calculated by reassembling the angular spectrum of all waves k after traversing the material. The focal volume is subsequently calculated by evaluating the diffraction integral in Eq. (1) with the modified wavefront.

For a fixed object expressed in angular frequency space, the EMC method allows a direct calculation of the distorted wavefront, and thus the simulation of the perturbed electric field in the vicinity of the focal volume. However, since different objects produce different wavefronts, the resulting focal volumes can vary broadly as a function of the shape, density, position, and refractive index properties of the scattering objects in the sample. The separate evaluation of each particular arrangement of scattering objects is not a convenient approach to distill the general trends of changes to the focal volume due to scattering. Nonetheless, the EMC method can be used in an approximate approach for determining the effects of scattering on the focal fields. To this end, the wavefront is determined by considering the effects of random scattering events on the amplitude and phase of the angular spectrum. Instead of evaluating fixed scattering objects, each wave k is randomly scattered in the medium through Monte Carlo sampling of uncorrelated scattering events. The effective electric field is then found by the coherent superposition of the scattered waves. Such a coherent summation approach has been used to calculate the effective field at selected points in the focal volume [35]. Here we use a similar calculation to determine the effective field of a given wave vector k in the angular spectrum representation. The amplitude and phase of the field associated with the plane wave exiting at a given angle θ_p is then calculated by taking the coherent sum of the contributions with the same exit angle:

Re (E_{∥ p}) = \frac{1}{\sqrt{N_{p} S_{p}}} Σ_{j = 1}^{N_{p}} Re (E_{∥ p, j})

Im (E_{∥ p}) = \frac{1}{\sqrt{N_{p} S_{p}}} Σ_{j = 1}^{N_{p}} Im (E_{∥ p, j})

where N_p is the number of waves detected in bin p, S_p is the surface area of bin p, and E _{∥p, j} is the jth wavefront detected in bin p whose amplitude is normalized by the total number of waves launched. The final amplitude is determined by $E_{∥ p} = \sqrt{Re {(E_{∥ p})}^{2} + Im {(E_{∥ p})}^{2}}$ . Similar expressions are used to calculate E _⊥p. Note that the coherent sum acts as a coherent filter, i.e., angular components that exhibit a large phase variation, due to random scrambling of the phase, have smaller final amplitudes than angular components where the phase is more conserved [35]. Thus we will refer to E _∥p as the coherent amplitude. Consequently, the random scattering approximation allows the calculation of an effective wavefront with an amplitude and phase resulting from the summation over many random scattering events. We use the procedure outlined above to extract general trends in the amplitude and phase of the focal fields in media with different scattering properties.

3. Methods

We examine tissue slabs with optical absorption and scattering coefficients of µ_a = 0.02/mm and µ′_s = 2/mm which are representative of human dermis at λ = 800 nm [44, 45]. Detailed scattering characteristics were determined using Mie theory [43] with spherical scatterers with a relative refractive index of 1.035. Spheres of radius a=0.2961 µm, 0.1873 µm, and 0.001 µm were used to produce anisotropy coefficients g = 0.8, 0.6, and 0, respectively. The number density of scatterers was adjusted to provide a transport mean free path of l* = 1/(µ′_s + µ_a) = 495 µm for all samples.

We modeled slab thicknesses of T = 0–1.5l*. The Monte Carlo simulations generated G(θ, ϕ → θ′,ϕ′) by launching 10⁹–2×10¹⁰ wavefronts for each slab thickness which produced relative errors of <0.1% in the wave count tallies. The focal volume was discretized into 61×61×121 voxels (x,y,z), over a focal volume that measures 3µm laterally and 6µm axially for capturing the focal field. The wave vector of each initial wavefront was selected from a Lambertian distribution sampled over 31 angular bins, while the exiting wavefront was tallied in 31θ bins in accordance with the solid angle of the objective lens. All simulations launched wavefronts of x-polarized light defined by E _∥ = 1 and E _⊥ = 0 into the slab. The EMC simulation provides the transfer function G(θ, ϕ → θ′,ϕ′) by providing the effective field for a particular exiting angle (θ′,ϕ′) given all incident angles (θ, ϕ). The focal volumes were constructed by computing Eqs. (1)–(3) numerically for numerical apertures of NA = 0.81, 1.16, and 1.31.

4. Results and Discussion

4.1. Spatial broadening of the focal volume

We first examine the effect of slab thickness on the spatial dispersion and strength of the focal field. We considered three slabs types all of which have a fixed transport mean free path l* = 495µm but have varying single-scattering anisotropy coefficients of g = 0,0.6,0.8. Holding l* constant for these three slab types resulted in µ_s values of 2.0, 5.0, 10.0mm⁻¹, respectively. In Fig. 2 we plot the variation of the (a) lateral and (b) axial dimensions of the focal field i.e., the measured full width at half maximum (FWHM), with the slab thickness for a numerical aperture (NA) of 1.16. These values are normalized to those predicted by diffraction alone in the absence of scattering.

These results demonstrate that for turbid media with equivalent values for l*, increases in g produce much stronger axial and lateral dispersion of the focal field. This occurs even for slab thicknesses exceeding 4l*; a thickness where one might expect diffusive light transport to be operative. However, because the focal field is formed primarily by the wavefront components that remain ‘in-phase’ after propagating through the material, the larger single-scattering coefficients µ_s associated with higher g result in a more pronounced spatial broadening of the focal volume. Although the single scattering phase function is more forward-directed for higher g, these plots suggest that the broadening of the focal fields may be governed predominantly by µ_s rather than the scattering direction. In addition, comparison of Fig. 2a with Fig. 2b reveals a stronger dispersion along the axial dimension. This is an expected result, as it is known from diffraction theory that the field confinement in the axial dimension is more sensitive to field aberrations than the field distribution in the lateral dimension [46, 47].

Fig. 2. (a) Lateral and (b) axial dimension of the focal field (FWHM) as a function of slab thickness T in units of l* for anisotropy coefficients g = 0, 0.6, 0.8 for a numerical aperture NA = 1.16. l* = 495µm in all samples.

Download Full Size | PDF

To illustrate the primacy of single-scattering in governing the spatial dispersion of the focal field, in Figs. 3(a) and (b), we plot the lateral and axial widths (FWHM) of the focal field versus the slab thickness T expressed in multiples of the single-scattering mean free path, l_s = 1/µ_s for three numerical apertures: 0.81, 1.16, and 1.31. When plotted in this fashion, the broadening characteristics of the focal field for a given numerical aperture in all three slabs fall onto a single curve. This confirms that the broadening of the focal field is governed solely by the expected number of scattering interactions independent of the angular distribution of the single-scattering phase function. Hence, for fixed values of l*, µ_a, and NA, the depth-dependent broadening of the focal field is given by a single ‘master curve’ when expressing the slab thickness in multiples of the single-scattering mean free path l_s.

Figures 3(a) and (b) also show the dependence of focal field broadening on the NA of the focusing lens. As one might expect, the relative broadening of the focal field is more significant for larger numerical apertures with increasing material thickness. This is because the formation of a diffraction-limited focal volume at higher numerical apertures relies on the unimpaired propagation of larger θ (off-axis) components of the wavefront which, due to their longer mean propagation distance through the slab, are more vulnerable to scattering or phase delay as they traverse the slab [25]. This observed NA dependence of focal volume dispersion is in line with experimental results [48] and incoherent MC calculations [27].

Note that although broadening of the focal fields is observed in these simulations, the overall effect of scattering on the confinement of the focal field distribution is modest. For instance, for a focusing depth corresponding to two single-scattering mean free paths (l_s = 2) and a numerical aperture of 1.16, the lateral broadening is less than 10%. This modest broadening agrees with reported multiphoton microscopy experiments that examine the quality of tightly focused excitation volumes in scattering media. These studies suggest that although the loss of the excitation amplitude deeper into the medium can be severe, the corresponding broadening of the focal excitation volume is generally modest [25, 49].

Fig. 3. (a) Lateral and (b) Axial FWHM as a function of slab thickness in terms of l_s for anisotropy coefficients g = 0,0.6,0.8 and numerical apertures NA = 0.81 (○), 1.16 (△), 1.31 (□).

Download Full Size | PDF

4.2. Loss of coherent amplitude

We next examine the attenuation of coherent amplitude with slab thickness. In Figs. 4(a) and (b) we plot the maximum amplitude in the focal volume as a function of slab thickness normalized with respect to the (a) transport mean free path l* and (b) single scattering mean free path l_s, respectively. Here the electric field amplitude is normalized relative to the maximum amplitude in the focal volume in the absence of scattering. Figure 4(a) demonstrates that for slabs of equivalent thickness (since l* = 495µm in all the slabs), the attenuation of the focal field is more severe for slabs with higher g which, accordingly, have a higher scattering coefficient µ_s.

Moreover, we find that when comparing amongst slabs with a fixed value of g, the attenuation becomes more significant for higher NAs. This latter observation can be rationalized by the more pronounced scattering-induced amplitude loss of field components traveling with larger θ when focusing with a higher NA lens. These features indicate, as was the case for the spatial dispersion of the focal field, that single scattering governs the attenuation of the focal field. Similarly, when plotting these results in terms of the single-scattering mean free path l_s in Fig. 4(b), the depth dependent field attenuation effectively falls onto a single curve with relatively small variations due to differences in g and NA.

Fig. 4. Normalized maximum amplitude as a function of slab thickness in terms of (a) l* and (b) l_s for numerical apertures NA = 0.81, 1.16, 1.31 and anisotropy coefficients g = 0(○), 0.6(△), 0.8(□).

Download Full Size | PDF

Recall that in the cases studied here, we compute the CADF by launching a Lambertian distribution of wavefronts with parallel polarization and capture the angular dispersion of each electric field component that is produced as a result of propagation through the slab. In Fig. 5 we plot the amplitude of the CADF for slab thicknesses of 1, 3, and 5 l_s for a single scattering anisotropy of g = 0.8. These plots clearly show increased attenuation and angular dispersion of the coherent amplitude with increasing angle of incidence and slab thickness.

Fig. 5. CADF for slab thicknesses of (a) 1l_s, (b) 3l_s, and (c) 5l_s for a single scattering anisotropy of g = 0.8. θ and θ′ represent the incident and exiting angle of the wavefront, respectively. The color bar represents a logarithmic scale.

Download Full Size | PDF

However, the formation of the focal field is governed solely by the amplitude, phase, and angular distribution of the wavefronts that exit the slab independent of the incident angles of the wavefront. In Fig. 6 we plot the normalized amplitude of the parallel component of the electric field E _∥ as a function of the exiting angle for slab thicknesses T = 1, 3, and 5l_s and fixed single-scattering anisotropy g = 0.8. These results are generated by integrating the CADF over θ. This plot reveals the overall attenuation of E _∥ with slab thickness. Moreover, this plot shows that the effect of increasing the normalized slab thickness is akin to a low-pass angular coherence filter. It is this low-pass filtering effect that is the origin of the spatial dispersion of the focal volume with depth. Hence, analysis of the CADF has allowed us to link the observed attenuation and spatial broadening of the focal fields to the angular distribution of the exiting wavefront as predicted by the EMC simulations.

Fig. 6. Normalized E _∥ amplitude versus exiting angle θ′ for slab thicknesses of T = 1l_s(○),3l_s(△),5l_s(□) for a single scattering anisotropy g = 0.8.

Download Full Size | PDF

4.3. Depolarization of the focal field

An important feature of this EMC approach is the ability to examine the effects of slab thickness and scatterer type on the depolarization that results from wavefront propagation to the focal volume. One motivation for this work is to provide a computational framework to predict the resolution and depth of linear and nonlinear optical microscopy techniques currently in use for imaging of thick tissues. A question relevant to polarization-sensitive methods such as harmonic generation microscopy and the co-registration of images relative to methods such as multiphoton microscopy is how the amplitude loss of incident field component E _∥ might relate to the appearance/generation of the E _⊥ components i.e., to the depolarization of the incident light.

To examine the interrelationship between loss of coherent amplitude and depolarization, in Fig. 7(a) we plot the normalized electric field component E _∥ as a function of the exit angle θ′ for slab thicknesses T = 1, 3, and 5l_s and single-scattering anisotropy g = 0, 0.6, and 0.8. As was the case in Figs. 3–5, we find a similar behavior for transmitted electric field in slabs with different g but identical thickness relative the the single-scattering length l_s. As in Fig. 6 we see the low pass filter effect of slabs with larger thickness. To compare, we plot in Fig. 7(b) a measure of the depolarization, D = 1−(E _⊥ ⁄ E _∥), as a function of θ′ in these same slabs. Figure 7(b) displays characteristics similar to Fig. 6 in that increasing slab thicknesses act as a low pass filter, this time for polarization as opposed to coherent amplitude. This agrees with experimental measurements of media with tissue-like properties which showed that the rate of depolarization of the incident radiation appears to be relatively small over length scales relevant to focused light [50]. However, comparison of the Figs. 7(a) with 7(b) reveals that propagation through a turbid slab of a fixed thickness consistently provides a more stringent filter for coherent amplitude as compared to polarization. This suggests that, for the scattering media examined, when using linearly-polarized incident light, the focal volume is formed predominantly by light with the same polarization state.

5. Conclusion

In this work we have used plane wave electric field Monte Carlo (EMC) simulations combined with an angular spectrum representation diffraction theory for focused fields to gain insight in the effects of tissue scattering on the formation of a tightly focused laser spot in turbid media. Unlike incoherent Monte Carlo models, the EMC approach preserves the wave properties of the focused radiation and is thus better suited to study the influence of scattering on the diffraction limited focal volume. We used this model in the random scattering limit to establish a fundamental link between the macroscopic tissue scattering parameters and the relative amplitude, spatial broadening, and depolarization of the focal field.

Fig. 7. (a) E _∥ normalized by maximum value, and (b) (1−E _⊥/E _∥), as a function of the wave’s exiting angle (θ′) for anisotropy coefficients g = 0,0.6,0.8 and slab thicknesses of T = 1l_s(○),3l_s(△),5l_s(□).

Download Full Size | PDF

Our results indicate that the loss of focal field amplitude with focusing depth is governed predominantly by the scattering coefficient µ_s rather than by the scattering anisotropy g. This implies that the properties of the focal volume are dictated primarily by the number of scattering events and that the direction of scattering plays a relatively minor role. The major effect of tissue scattering is the increased amplitude loss and phase delay for wavefront components k with large propagation angles. This modifies the angular spectrum, in that the low angle k components gain more importance relative to the high angle components. Effectively, the scattering medium acts as a low pass filter of the angular spectrum, and results in a broadening of the focal fields in both the lateral and axial dimensions. This low pass filtering mechanism also explains why effective broadening of the focal fields is less severe for lower NA objective lenses, as the angular spectrum of a low NA lens intrinsically encompasses only lower angular components. Although effective broadening of the focal volume results directly from tissue scattering, our model confirms quantitatively that spatial distortions of the focal fields are relatively minor.

In addition, the full vectorial nature of the EMC approach allows a direct assessment of the depolarization rate of the focal fields. Our simulations indicate that while the depolarization of higher angular components can be substantial for large slab thicknesses l*, this effect is subordinate to the corresponding coherent amplitude loss. We conclude that for the scattering media examined, the loss of coherent amplitude due to phase scrambling is much more severe than depolarization within the focal field of a tightly focused laser beam in turbid media.

All of the trends predicted by the simulations presented here are in good agreement with experimental observations. The EMC model for focused light in the limit of the random scattering approximation thus provides a reliable prediction of the excitation field in turbid media based on macroscopic tissue scattering parameters. We expect this approach to be particularly relevant to predict the performance of coherent imaging techniques applied to turbid media, which crucially relies on a full assessment of the amplitude and phase of the focal volume.

Acknowledgements

We acknowledge support from the Laser Microbeam and Medical Program (LAMMP) a NIH Biomedical Technology Resource Center (P41-RR01192). CKH acknowledges support from the National Institutes of Health (NIH, K25-EB007309) and EOP acknowledges support from the National Science Foundation (NSF, CHE-0847097).

References and links

1. V. Tuchin, Tissue Optics (SPIE Press, 2007).

2. P. Theer, M. T. Hasan, and W. Denk, “Two-photon imaging to a depth of 1000 µm in living brains by use of a Ti:Al₂O₃ regenerative amplifier,” Opt. Lett. 28(12), 1022–1024 (2003). [CrossRef] [PubMed]

3. M. Balu, T. Baldacchini, J. Carter, T. B. Krasieva, R. Zadoyan, and B. J. Tromberg, “Effect of excitation wavelength on penetration depth in nonlinear optical microscopy of turbid media,” J. Biomed. Opt. 14(1), 010508 (2009). [CrossRef] [PubMed]

4. M. J. Booth, M. A. A. Neil, R. Juskaitis, and T. Wilson, “Adaptive aberration correction in a confocal microscope,” Proc. Natl. Acad. Sci. U.S.A. 99(9), 5788–5792 (2002). [CrossRef] [PubMed]

5. L. Sherman, J. Y. Ye, O. Albert, and T. B. Norris, “Adaptive correction of depth-induced aberrations in multiphoton scanning microscopy using a deformable mirror,” J. Microsc. 206(1), 65–71 (2002). [CrossRef] [PubMed]

6. T. M. Nieuwenhuizen, A. Lagendijk, and B. A. van Tiggelen, “Resonant point scatterers in multiple scattering of classical waves,” Phys. Lett. A 169(3), 191–194 (1992). [CrossRef]

7. A. K. Dunn and R. Richards-Kortum, “Three-dimensional computation of light scattering from cells,” IEEE J. Sel. Top. Quantum Electron. 2(4), 898–905 (1996). [CrossRef]

8. R. Drezek, A. Dunn, and R. Richards-Kortum, “Light scattering from cells: finite-difference time-domain simulations and goniometric measurements,” Appl. Opt. 38(16), 3651–3661 (1999). [CrossRef] [PubMed]

9. C. Liu, C. Capjack, and W. Rozmus, “3-D simulation of light scattering from biological cells and cell differentiation,” J. Biomed. Opt. 10(1), 014007 (2005). [CrossRef] [PubMed]

10. I. R. Çapoglu, A. Taflove, and V. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express 16(23), 19208–19220 (2008). [CrossRef] [PubMed]

11. M. S. Starosta and A. K. Dunn, “Three-dimensional computation of focused beam propagation through multiple biological cells,” Opt. Express 17(15), 12455–12469 (2009). [CrossRef] [PubMed]

12. A. Ishimaru, Wave Propagation and Scattering in Random Media, Vols. I and II (Academic Press, 1978).

13. L. Tsang, J. A. Kong, and K. H. Ding, Scattering of Electromagnetic Waves: Theories and Applications (Wiley, 2000).

14. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20(1), 92–98 (2003). [CrossRef] [PubMed]

15. G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. 7(8), 1519–1527 (1968). [CrossRef] [PubMed]

16. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). [CrossRef] [PubMed]

17. I. Lux, and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (CRC Press, 1991).

18. X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt. 7(3), 279–290 (2002). [CrossRef] [PubMed]

19. J. S. You, C. K. Hayakawa, and V. Venugopalan, “Frequency domain photon migration in the δ- P1 approximation: analysis of ballistic, transport, and diffuse regimes,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72(2), 021903 (2005). [CrossRef] [PubMed]

20. J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A 13(5), 952–961 (1996). [CrossRef]

21. C. M. Blanca and C. Saloma, “Monte carlo analysis of two-photon fluorescence imaging through a scattering medium,” Appl. Opt. 37(34), 8092–8102 (1998). [CrossRef] [PubMed]

22. 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(13), 2944–2949 (1999). [CrossRef] [PubMed]

23. L. V. Wang and G. Liang, “Absorption distribution of an optical beam focused into a turbid medium,” Appl. Opt. 38(22), 4951–4958 (1999). [CrossRef] [PubMed]

24. X. S. Gan and M. Gu, “Effective point-spread function for fast image modeling and processing in microscopic imaging through turbid media,” Opt. Lett. 24(11), 741–743 (1999). [CrossRef] [PubMed]

25. A. K. Dunn, V. P. Wallace, M. Coleno, M. W. Berns, and B. J. Tromberg, “Influence of optical properties on two-photon fluorescence imaging in turbid samples,” Appl. Opt. 39(7), 1194–1201 (2000). [CrossRef] [PubMed]

26. X. Deng, X. Gan, and M. Gu, “Monte Carlo simulation of multiphoton fluorescence microscopic imaging through inhomogeneous tissuelike turbid media,” J. Biomed. Opt. 8(3), 440–449 (2003). [CrossRef] [PubMed]

27. 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(16), 3321–3329 (2003). [CrossRef] [PubMed]

28. 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(2), 024013 (2006). [CrossRef] [PubMed]

29. 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(1), 269–278 (2007). [CrossRef]

30. J. M. Schmitt and A. Knuttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14(6), 1231–1242 (1997). [CrossRef]

31. D. J. Smithies, T. Lindmo, Z. Chen, J. S. Nelson, and T. E. Milner, “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol. 43(10), 3025–3044 (1998). [CrossRef] [PubMed]

32. A. Tycho, T. M. Jørgensen, H. T. Yura, and P. E. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt. 41(31), 6676–6691 (2002). [CrossRef] [PubMed]

33. G. Xiong, P. Xue, J. Wu, Q. Miao, R. Wang, and L. Ji, “Particle-fixed Monte Carlo model for optical coherence tomography,” Opt. Express 13(6), 2182–2195 (2005). [CrossRef] [PubMed]

34. D. G. Fischer, S. A. Prahl, and D. D. Duncan, “Monte Carlo modeling of spatial coherence: free space diffraction,” J. Opt. Soc. Am. A 25(10), 2571–2581 (2008). [CrossRef]

35. V. R. Daria, C. Saloma, and S. Kawata, “Excitation with a focused, pulsed optical beam in scattering media: diffraction effects,” Appl. Opt. 39(28), 5244–5255 (2000). [CrossRef] [PubMed]

36. M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express 12(26), 6530–6539 (2004). [CrossRef] [PubMed]

37. K. Phillips, M. Xu, S. K. Gayen, and R. R. Alfano, “Time-resolved ring structure of circularly polarized beams backscattered from forward scattering media,” Opt. Express 13(20), 7954–7969 (2005). [CrossRef] [PubMed]

38. J. Sawicki, N. Kastor, and M. Xu, “Electric field Monte Carlo simulation of coherent backscattering of polarized light by a turbid medium containing Mie scatterers,” Opt. Express 16(8), 5728–5738 (2008). [CrossRef] [PubMed]

39. C. K. Hayakawa, V. Venugopalan, V. V. Krishnamachari, and E. O. Potma, “Amplitude and phase of tightly focused laser beams in turbid media,” Phys. Rev. Lett. 103(4), 043903 (2009). [CrossRef] [PubMed]

40. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems 2: structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A 253(1274), 358–379 (1959). [CrossRef]

41. L. Novotny, and B. Hecht, Principles of Nano-Optics (Cambridge University Press, 2006).

42. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001). [CrossRef] [PubMed]

43. C. F. Bohren, and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley and Sons, 1983).

44. T. L. Troy and S. N. Thennadil, “Optical properties of human skin in the near infrared wavelength range of 1000 to 2200 nm,” J. Biomed. Opt. 6(2), 167–176 (2001). [CrossRef] [PubMed]

45. S. L. Jacques, “Skin Optics,” http://omlc.ogi.edu/news/jan98/skinoptics.html.

46. B. R. A. Nijboer, “The diffraction theory of optical aberrations. Part I: General discussion of the geometrical aberrations,” Physica 10(8), 679–692 (1943). [CrossRef]

47. T. Wilson and A. R. Carlini, “Aberrations in confocal imaging systems,” J. Microsc. 154, 243–256 (1998).

48. C. K. Tung, Y. Sun, W. Lo, S. J. Lin, S. H. Jee, and C. Y. Dong, “Effects of objective numerical apertures on achievable imaging depths in multiphoton microscopy,” Microsc. Res. Tech. 65(6), 308–314 (2004). [CrossRef] [PubMed]

49. C. Y. Dong, K. Koenig, and P. So, “Characterizing point spread functions of two-photon fluorescence microscopy in turbid medium,” J. Biomed. Opt. 8(3), 450–459 (2003). [CrossRef] [PubMed]

50. N. Ghosh, H. S. Patel, and P. K. Gupta, “Depolarization of light in tissue phantoms - effect of a distribution in the size of scatterers,” Opt. Express 11(18), 2198–2205 (2003). [CrossRef] [PubMed]

Electric field Monte Carlo simulations of focal field distributions produced by tightly focused laser beams in tissues

Abstract

1. Introduction

2. Theory

2.1. Angular spectrum representation of diffraction theory

2.2. Monte Carlo simulations

3. Methods

4. Results and Discussion

4.1. Spatial broadening of the focal volume

4.2. Loss of coherent amplitude

4.3. Depolarization of the focal field

5. Conclusion

Acknowledgements

References and links

Cited By

Figures (7)

Equations (17)

Biomedical Optics Express