Modeling the tight focusing of beams in absorbing media with Monte Carlo simulations

Arnd R. Brandes; Ahmed Elmaklizi; H. Günhan Akarcay; Alwin Kienle

doi:10.1117/1.JBO.19.11.115003

13 November 2014 Modeling the tight focusing of beams in absorbing media with Monte Carlo simulations

Arnd R. Brandes, Ahmed Elmaklizi, H. Günhan Akarcay, Alwin Kienle

Author Affiliations +

Funded by: Evangelisches Studienwerk e.V. Villigst, the Baden-Württemberg Stiftung and the Swiss National Science Foundation

Journal of Biomedical Optics, Vol. 19, Issue 11, 115003 (November 2014). https://doi.org/10.1117/1.JBO.19.11.115003

Abstract

A severe drawback to the scalar Monte Carlo (MC) method is the difficulty of introducing diffraction when simulating light propagation. This hinders, for instance, the accurate modeling of beams focused through microscope objectives, where the diffraction patterns in the focal plane are of great importance in various applications. Here, we propose to overcome this issue by means of a direct extinction method. In the MC simulations, the photon paths’ initial positions are sampled from probability distributions which are calculated with a modified angular spectrum of the plane waves technique. We restricted our study to the two-dimensional case, and investigated the feasibility of our approach for absorbing yet nonscattering materials. We simulated the focusing of collimated beams with uniform profiles through microscope objectives. Our results were compared with those yielded by independent simulations using the finite-difference time-domain method. Very good agreement was achieved between the results of both methods, not only for the power distributions around the focal region including diffraction patterns, but also for the distribution of the energy flow (Poynting vector).

1. Introduction

The generation of localized light fields through high numerical aperture ( $NA \geq 0.7$ ¹) optical systems has become increasingly common in the field of biophotonics, with applications ranging from optical tweezers to microscopy and spectroscopy.²^–⁶ Yet, such systems ask for reliable simulation tools, not only for the design of the optics, but also for a better understanding of light-matter interactions (e.g., absorption and scattering) within the focal volume. The tight focusing of beams induces diffraction patterns in the focal volume where the polarization of the illumination beam plays an important role.⁷^,⁸ Consequently, if the tight focusing of beams is to be modeled, the consideration of the wave nature of the light becomes paramount. This can be successfully achieved by solving numerically Maxwell equations, e.g., by carrying out finite-difference time-domain (FDTD) simulations⁹ for small volumes.

The Monte Carlo (MC) algorithm¹⁰ is commonly used for calculations of the light propagation in scattering and absorbing media, and the problem of focused beams has been assessed by various groups and is still a matter of ongoing research.¹¹^–¹⁴ We present here an alternative MC approach that enables the handling of the diffraction within absorbing media introduced by a coherent light source. Moreover, the high sampling number required for a destructive interference of a coherent source in a vectorial MC¹⁵ is rendered obsolete by this approach. In this contribution, we investigate the feasibility of including diffraction of the incident light into MC simulations in order to model tightly focused beams and Fraunhofer diffraction patterns in arbitrarily large simulation volumes in a relatively simple manner.

This paper is organized as follows: in Sec. 2, we describe how the MC simulations are combined with the so-called angular spectrum of the plane waves technique (ASPWs) to simulate the light propagation and to introduce, at the same time, our “direct extinction method.” In Sec. 3, we present the results yielded by the developed method and compare them with those obtained from FDTD simulations, both in vacuo and in absorbing media. We are interested not only in the power distribution over the sample, but also in the distribution of the electric and magnetic field components and, subsequently, in the energy flow (represented by the Poynting vector) at every point. In the last section, we discuss the obtained results and present conclusions on the feasibility of our approach.

2. Theory

We consider the tight focusing of a monochromatic beam through a microscope objective (see Fig. 1). A collimated laser beam with a uniform profile and propagating along the $x$ -direction is focused by an aplanatic cylindrical lens.¹⁶ The focal volume lies within the simulation grid. To simplify the computational effort and the simulation of larger grids, we confine our investigation here to the two-dimensional (2-D) case, so that only the $x, y$ -plane is considered. We illustrate our purpose with two types of polarizations for the illuminating beam: first, linearly polarized along the $z$ -axis, and second, radially polarized (within the $x, y$ -plane). The sample and its surrounding medium are assumed to have the same index of refraction, i.e., there is no Fresnel reflection at the boundaries.

Fig. 1

Geometry of the simulations, as carried out in this paper. The setup is infinite in $z$ -direction.

We distinguish two cases: first, the beam propagates in vacuo, and second, the beam propagates in an absorbing, yet nonscattering material.

2.1.

Direct Extinction Method within the Monte Carlo Algorithm

When dealing with light propagation problems, the MC method constitutes a well-established numerical alternative to the analytical solutions of the radiative transfer equation. ¹⁰^,¹⁷^–¹⁹ The core idea consists of stochastically modeling the light propagation with individual photon paths akin to random walks. A single photon path is composed of successive ray segments, each of which is defined by (1) a starting point, (2) a propagation direction, (3) and a segment length. Since abundant documentation can be found on the matter,¹⁰^,¹⁵^,²⁰ we do not repeat the successive steps needed to generate such paths; instead, we emphasize the novelty of our approach. Traditionally, the initial starting point of the very first ray segment is chosen at the position of the illuminating source, or on the surface of the sample, from where the illumination is performed. Here, we propose an alternative method where the initial starting point of the rays coincides to a “target position”. This position is independently sampled from a probability distribution calculated with the ASPW method and which corresponds to the probability that the ray has not interacted (by scattering or absorption) with the sample yet. We confine this study to absorbing materials as a preliminary investigation.

In this section, we detail the calculation of the probability distribution for the incident beam and the inclusion of the absorption into the ASPW method. Further, we explain how our method can be used to calculate the local Poynting vectors, as these can help to model scattering processes at a later stage. In all of these simulations, only the 2-D case is studied to simplify the computation and to decouple the polarization of the light into two orthogonal components: parallel and perpendicular. However, the consideration of the 2-D case is beyond a mere simplification of the three-dimensional (3-D) case. “In fact, both the analyses and the final results are sufficiently different in the two cases to call for separate treatments.”²¹

2.2.

Angular Spectrum of Plane Waves Method in Vacuo

In order to calculate the power distribution corresponding to the probability of finding a photon propagating in the simulation space, we need to know the electric field $E (r)$ at every point $r = (x, y)$ . This can be obtained with the ASPW method, which is particularly well suited for problems involving arbitrarily shaped beams where the paraxial approximation fails, as in our current proposal.¹⁶^,²² The ASPW method is a mathematical technique that allows for the representation of complex beam profiles by a coherent superposition of plane waves: this is achieved knowing that the field distribution in the focal plane is given by the Fourier transform of the field distribution in the far field, i.e., at the lens’ surface (we neglect the contribution of the evanescent waves; as in our case, the focal plane is assumed to be located sufficiently afar from the lens).

Since, in a first step, the monochromatic beam is propagating in the $x$ -direction in a nonabsorbing and nonscattering 2-D space, each plane wave in the $x, y$ -plane is characterized by a wave vector $k$ , which is defined by the wave’s incidence angle $θ$ onto the focal plane (see Fig. 1). Thus, if we consider a linearly polarized beam such that its electric field is oscillating along the $z$ -direction ( $E = E_{z}$ ), the angular spectrum representation of the beam yields the expression of the field

Eq. (1)

E_{z} (r) = \frac{1}{2 π} \int_{- \infty}^{\infty} E_{z}^{0} (k_{y}) \exp [j (x - x_{ph}) \sqrt{k^{2} - k_{y}^{2}} + j k_{y} y] d k_{y},

at every point

r

in the simulation space where

k = \sqrt{k_{x}^{2} + k_{y}^{2}}

is the wave number in the medium and

x_{ph}

is the position of the focal plane (see Fig. 2). As for

E_{z}^{0} (k_{y})

, it expresses the field in the focal plane and is obtained by the transformation of a collimated beam with a uniform profile truncated by an aperture incident on the cylindrical lens. We assume the field to be zero outside the aperture (Kirchhoff approximation). Further, we follow here the reasoning of Richards, Wolf, and Born,¹⁶^,²³ where the intensity law of geometrical optics allows for the calculation of the fields after passing through the lens, which leads to the formulation²⁴ of

Eq. (2)

E (r) = \frac{k}{2 π} \int_{- θ_{\max}}^{θ_{\max}} \hat{E} a (r, θ) d θ,

where

\hat{E}

denotes the amplitude of the incident electric field. Each plane wave is scaled by its respective factor of

a (r, θ)

to give a linearly polarized beam defined as

Eq. (3)

a (r, θ) = (\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}) \cdot \exp [j k ((x - x_{p h}) \cos θ + y \sin θ)] \sqrt{\cos (θ)} .

Fig. 2

Comparison between the grid sizes of FDTD and MC.

Note that the Kirchhoff approximation, combined with the assumption that the focal region is located sufficiently afar from the lens, corresponds to the so-called Debye-approximation.²¹^,²⁵ In order to implement a radially polarized beam, $a (r, θ)$ needs to be modified:²⁶

Eq. (4)

a_{rad} (r, θ) = (\begin{matrix} \sin θ \\ \cos θ \\ 0 \end{matrix}) sgn (θ) a_{z} (r, θ) .

The power $P_{v} (r) \propto E (r) \cdot E^{*} (r)$ leads to a probability density function $p_{v} (r)$ , which expresses the probability of finding a photon at a position $r$

Eq. (5)

p_{v} (r) = \frac{P_{v} (r)}{\int_{r} P_{v} (r) d r} .

2.3.

Modeling the Absorption

Now consider the beam propagating in an absorbing material with a refractive index $n = n + j κ$ , such that its absorption coefficient is defined by $μ_{a} = 2 k κ$ .

If we denote the fields of the beam in vacuo by $E (r)$ , the corresponding expression of the fields in the absorbing material is $E^{e} (r, μ_{a})$ , where the superscript “e” refers to the extinction. The influence of the medium can be represented by inhomogeneous plane waves²³ composing the beam, based on Beer-Lambert law’s exponential decay (we neglect the Fresnel reflection which occurs between the nonabsorbing and the absorbing medium). The modified ASPW representation is defined as

Eq. (6)

E^{e} (r, μ_{a}) = \int_{- θ_{\max}}^{θ_{\max}} \hat{E} a^{e} (r, θ, μ_{a}) d θ,

where the integrand introduces the Beer-Lambert law as

Eq. (7)

a^{e} (r, θ, μ_{a}) = a (r, θ) \exp [- \frac{1}{2} μ_{a} \frac{x}{\cos (θ)}] .

Analogously to the case in vacuo [Eq. (5)], the power $P_{a} (r) \propto E^{e} (r, μ_{a}) \cdot E^{e *} (r, μ_{a})$ leads to a probability density function $p_{a} (r)$ , which expresses the the probability of finding a photon at a position $r$ before it gets absorbed in the sample:

Eq. (8)

p_{a} (r) = \frac{P_{a} (r)}{\int_{r} P_{a} (r) d r} .

This constitutes the essence of the “direct extinction method,” which can later be extended to include scattering processes. The main advantage of this approach is that the diffraction phenomena can be recovered from the sampled probability distribution without having to implement the interference between photon paths.

For simplicity, in the following part of our manuscript, we shall use the term “beam” to refer to the beam which fulfills the assumptions described in this section and whose electric field is given by $E (r)$ in vacuo and by $E^{e} (r, μ_{a})$ in an absorbing medium.

2.4.

Ray Direction: Poynting Vector

As already mentioned, the present study is confined to nonscattering materials, i.e., the MC simulations do not need to propagate the photon paths, and the sampling from the aforesaid probability distribution suffices. However, if scattering processes are to be modeled, another key information besides the “first interaction point” needs to be known: the first path segment’s propagation direction.

To this end, we may utilize the energy flow or the energy flux density given by the local time-averaged Poynting vector $S (r)$

Eq. (9)

S (r) \propto real {E^{e} (r) \times H^{e *} (r)},

where the magnetic field

H^{e}

at a point

r

is calculated similarly to the electric field with the modified ASPW method to give

Eq. (10)

H^{e} (r, μ_{a}) = \int_{- θ_{\max}}^{θ_{\max}} b^{e} (r, θ, μ_{a}) d θ,

where

b

denotes the contribution of the single plane waves. For the linearly polarized beam we obtain

Eq. (11)

b^{e} (r, θ) = \hat{H} (\begin{matrix} \sin θ \\ \cos θ \\ 0 \end{matrix}) a_{z}^{e} (r, θ) .

$\hat{H}$ denotes the amplitude of the incident magnetic field. As for the radially polarized beam we obtain

Eq. (12)

b_{rad}^{e} (r, θ) = \hat{H} sgn (θ) a^{e} (r, θ) .

Note that for $μ_{a} = 0$ , the magnetic and electric fields $E^{e}$ and $H^{e}$ are reduced to the fields $E$ and $H$ in vacuo, respectively.

3. Simulations and Results

3.1.

Principle: in Vacuo

We simulated the 2-D version of the situation sketched in Fig. 1: the beam (as described in the previous section) propagating in the $x$ -direction is focused by an aplanatic microscope objective, so that the divergence angle is $θ_{\max} = 65 \deg$ , which lies beyond the paraxial limit. This corresponds to a numerical aperture of $NA = 0.9$ in air or $NA = 1.2$ in water, and is relatively common for modern confocal microscope objectives. As an example, we consider the monochromatic beam with a wavelength of $λ = 1 μ m$ . Note that by scaling, these results are also valid for any other wavelength.

The MC simulations were carried out in a square-shaped sample, with 200 bins in the $x$ - and $y$ -directions, with bin sizes of $Δ x = Δ y = λ / 25$ . As in the ASPW model, Fresnel reflection at the sample interfaces was neglected, i.e., the sample’s refractive index was assumed to be equal to that of the surrounding medium ( $n_{sample} = n_{air} = 1$ ). Since we are dealing with nonscattering materials, the MC simulations are confined to sampling “target positions” from the normalized probability distribution $p_{v} (r)$ calculated with the ASPW method, where 41 plane waves were superimposed. In vacuo, $p_{v} (r)$ simply expresses the probability of finding a photon at a position $r$ . The sampling is achieved with the rejection method. First, a pseudorandom position $r_{i}$ is drawn. (Pseudorandom numbers were generated with the Mersenne Twister²⁷ algorithm.) The corresponding probability density $p_{v} (r_{i})$ at that position is compared with a pseudorandom number $ξ \in (0,1]$ : If the number is such that $ξ > p_{v} (r_{i})$ , the position is rejected and the sampling is repeated.

In order to assess the accuracy of our results, we conducted independent FDTD simulations (more details on these simulations can be found elsewhere⁹^,²⁴^,²⁸). These simulations numerically solve the discretized time-dependent Maxwell’s equations of electrodynamics to propagate the electric and magnetic fields within a grid using a leap frog algorithm. As input, we need to specify for each plane wave its (1) incidence angle $θ$ on the sample, together with its subsequent starting position $r_{0} = (x = 0, y)$ , and (2) its electric field at $r_{0}$ . To remain consistent with the MC simulations, we propagated 41 plane waves with incidence angles lying in the range $[- θ_{\max}, + θ_{\max}]$ . For a given incidence $θ$ , an expression of the electric field at $r_{0}$ is obtained from the ASPW calculations [i.e., from Eq. (2)]. All FDTD simulations throughout the work were performed with a resolution of $λ / 25$ .

Results obtained with MC and FDTD simulations for a linearly polarized beam propagating in vacuo are shown in Fig. 3. In this paper, all simulation results are displayed after normalization by the probability density value at the focal point ( $x = x_{ph}$ , $y = 0$ ) of the simulation in vacuo, denoted by $p_{v, 0}$ and $p_{v, 0}^{'}$ for the linearly and radially polarized beams, respectively. Both types of simulations were performed with the same grid dimensions. The agreement between the two sets of data serves as a verification to the implementation of both the MC and the FDTD algorithms.

Fig. 3

Comparison of the probability distributions yielded by the Monte Carlo (MC) simulations (the number of sampled “target positions” are $10^{6}$ and $10^{9}$ on (a) and (b), respectively) and the finite-difference time-domain (FDTD) simulations for the linearly polarized beam focused in vacuo in the $x$ -direction. These distributions are displayed here in logarithmic scale over the range $8 λ \times 8 λ$ . The MC simulations are shown for $y / λ > 0$ and FDTD simulations for $y / λ \leq 0$ for both figures. The curves are cross sections of the probability profiles at the focal plane $x = x_{ph} = 4 λ$ .

Besides, as the FDTD method accounts for the near-field interactions, we can verify here the formation of the Fraunhofer diffraction patterns in the vicinity of the focal plane. However, we stress here the main constraint inherent to the MC technique: the number of samplings should be large enough to provide reasonable convergence.

3.2.

In Absorbing Media

As a following step, we compare the MC and FDTD simulations in the presence of absorption. This would be of interest, for instance, to investigate the excitation of fluorophores with confocal fluorescence systems. Again, the MC simulations were performed with a square-shaped sample of dimensions $8 λ \times 8 λ$ , such that the focal plane of the focused beam is located at $x_{ph} = 4 λ$ .

Four types of samples were considered, each sample having a different absorption coefficient ${(16 λ)}^{- 1} \leq μ_{a} \leq {(2 λ)}^{- 1}$ . Note that despite the large absorption values, we omit here re-emission, heat transfer, denaturation, etc., since these do not contribute to our present study. As opposed to the choice of these high absorption coefficients, it is also possible to simulate smaller values, but that requires a much larger grid in order to show an effect on the focus. However, an increased grid size is disadvantageous for the performance of the FDTD simulation.

We proceed according to our extinction method, by sampling “target positions” from $p_{a} (r)$ instead of $p_{v} (r)$ . These “target positions” bear a slightly different significance when absorption is involved: $p_{a} (r)$ represents the probability to find a photon in the sample before it gets absorbed along its path.

In the FDTD simulations, the absorption in the sample is taken into account by increasing the imaginary part of the refractive index: ${(64 π)}^{- 1} \leq κ \leq {(8 π)}^{- 1}$ . The grid spacing is preserved, but the overall dimensions of the simulation grid need to be modified. This is shown in Fig. 2. Precaution is needed when comparing the MC simulation results (based on a propagation space infinite in the $y$ -direction) with those yielded by FDTD simulations: in the case of larger divergence angles $θ_{\max}$ , the inclusion of absorption in the medium where the beam propagates requires an adjustment of the finite FDTD simulation grid. In the figure, the darker region corresponds to the MC simulation area. The symbol $x_{ph}$ indicates the position of the focal plane. The FDTD grid needs to be wide enough so that all plane waves start getting attenuated at $x = 0$ ; this principally concerns plane waves incident on the focal plane with a highly oblique angle $θ$ . The dashed area corresponds to the extra computational space: the choice of the additional space naturally affects the outcome of the simulations. The interaction between the incident plane waves (required by ASPW) and the absorbing medium requires an infinite boundary. However, FDTD can only calculate finite volumes.

An optimal grid size can be evaluated as shown in Fig. 4. Here, a larger grid space in the FDTD simulations compared with the grid considered in the MC simulations is evaluated. We show here the power distribution obtained with three different FDTD grid sizes for the linearly polarized beam focused in an absorbing material with $μ_{a} = {(2 λ)}^{- 1}$ in the $x$ -direction. The grid size in the FDTD simulations needs to be large enough to provide correct results by accounting for the absorption of plane waves with a wide angular incidence on the focal plane. The figure in the center and on the right reveals that the grid size needs to be above a certain threshold for the results to converge.²⁹

Fig. 4

Three different FDTD grid sizes displayed in logarithmic scale over the range $8 λ \times 8 λ$ for the linearly polarized beam, from (a) $x = 8 λ$ and $y = 8 λ$ , (b) $8 λ \times 16 λ$ , and (c) $8 λ \times 32 λ$ .

Simulation results obtained with the four types of samples are shown in Fig. 5 and more detailed profiles from the same results are shown in Fig. 6. It can be seen that good agreement is achieved between the MC and the FDTD simulations. Nonetheless, closer inspection reveals the appearance of discrepancies, both close to the sample’s illumination face and in the focal region (the latter case is not shown here but was observed, for instance, in the case of beams focused with a divergence angle of $θ_{\max} = 25 \deg$ ). Since similar discrepancies have been observed no matter what divergence angle $θ_{\max}$ was taken for the beam, the most plausible explanation for the differences is the presence of Fresnel reflection in the FDTD simulations. In the FDTD model, the plane waves propagate from air into a sample, whose refractive index has a nonzero imaginary component. However, as mentioned earlier, the MC simulations were implemented assuming that there is no refractive index mismatch along the plane waves’ path.

Fig. 5

Comparison of the probability distributions displayed here in logarithmic scale over the range $8 λ \times 8 λ$ and yielded by the MC simulations (with $10^{9}$ “target positions” sampled) and the FDTD simulations (grid size $8 λ \times 32 λ$ ) for the linearly polarized beam focused in the $x$ -direction in increasingly absorbing media. From (a) to (d): $μ_{a} = {(16 λ)}^{- 1}$ , ${(8 λ)}^{- 1}$ , ${(4 λ)}^{- 1}$ and ${(2 λ)}^{- 1}$ (the corresponding values for the imaginary component $κ$ of the material’s refractive index are ${(64 π)}^{- 1}$ , ${(32 π)}^{- 1}$ , ${(16 π)}^{- 1}$ , and ${(8 π)}^{- 1}$ , respectively). The curves are cross sections of the probability profiles at the focal plane $x = x_{ph}$ (MC in blue and FDTD in green).

Fig. 6

Quantitative representation of the results shown in Fig. 5. The probability profiles corresponding to the different absorbing materials are cross sections close to the illumination face (a), at $x = 0.01 λ$ and in the focal plane at $x = x_{ph} = 4 λ$ (b).

3.3.

Direction of the Poynting Vectors

In the last step of our investigation, we compared the orientation of the Poynting vectors $S (r)$ as calculated with both MC and FDTD simulations. In the MC method, the Poynting vector is calculated according to Sec. 2.4, averaged over a bin of subwavelength size $Δ x \times Δ y$ . On the other hand, in the FDTD simulations, the Poynting vector is determined at a certain grid point $(x, y)$ from the resulting electric and magnetic fields at the end of the simulation, see Eq. (9).

The results obtained by focusing a linearly polarized beam are shown in Fig. 7. Again, very good quantitative agreement is achieved between the two datasets. Yet, differences can be seen in the vicinity of the focal plane, in particular in the darker regions of the probability distribution where destructive interferences occur. In those regions, the orientation of the Poynting vector might vary considerably over a bin, in which case the averaging with the MC method does not yield identical results to the FDTD simulations. In other regions where discrepancies are not visible, the orientation of the Poynting vector remains relatively constant over a bin, and thus the averaging produces identical results to the FDTD simulations. Moreover, the effect of the aforementioned Fresnel reflection included in the FDTD calculations and neglected in the MC simulations should not be excluded in the interpretation of the differences. Avoiding the mismatch of the imaginary part of the refractive index in FDTD by adjusting the refractive index in the outer regions of the grid is a nontrivial task and would require further development of the FDTD algorithm, which is beyond the scope of the present feasibility study.

Fig. 7

Direction of the normalized Poynting vectors $S (r)$ plotted over probability distributions displayed here in logarithmic scale over the range $3 λ \times 3 λ$ , and calculated with both MC simulations ( $10^{9}$ “target positions” sampled) and FDTD simulations (grid size of $8 λ \times 32 λ$ ). (a) The linearly polarized beam is focused in vacuo. (b) The same beam is focused in an absorbing medium, where $μ_{a} = {(2 λ)}^{- 1}$ . In MC simulations, the binning size was kept at $Δ x = Δ y = λ / 25$ .

Additionally, results yielded by the focusing of a radially polarized beam are shown in Fig. 8: The corresponding Fraunhofer diffraction pattern noticeably differs from that of the previous example. As for the Poynting vectors in Fig. 8, the remarks on Fig. 7 still hold true. Here the observed spacial distribution depends on the resulting electric field from both $E_{x}$ and $E_{y}$ . Furthermore, the tight focusing of the radially polarized beam leads to a strong longitudinal electric field component in the $x$ -direction in the focal region, and the Poynting vector vanishes [as can be seen from the FDTD simulations of Fig. 8(a)].³⁰^–³²

Fig. 8

Radially polarized beam is focused in vacuo.(a) Distribution of the non-normalized Poynting vectors $S (r)$ plotted over probability distribution displayed here in logarithmic scale over the range $4 λ \times 4 λ$ , and calculated with FDTD simulations (grid size of $8 λ \times 32 λ$ ). (b) Same distribution of the normalized Poynting vectors obtained with both MC simulations ( $10^{8}$ “target positions” sampled) and the FDTD simulations, plotted over probability distribution displayed here in logarithmic scale over the range $8 λ \times 8 λ$ .

4. Discussion and Conclusions

In this contribution, we have investigated the feasibility of modeling the tight focusing of beams with MC simulations based on the “direct extinction method,” where “target positions” for the photon paths are sampled from probability distributions separately calculated with the ASPW technique. The simulations present an efficient tool to examine some of the wave properties of the light (Fraunhofer diffraction pattern, polarization state, and the direction of the energy flow) in the focal region of a high-aperture system, and unlike the FDTD simulations, the MC method presented here does not impose limitations on the dimensions of the observed sample.

Naturally, the accuracy of the results depends not only on the spatial resolution, but also on the number of sampled “target positions” and the number of plane waves used. While the FDTD grid has to be fine enough to avoid numerical artifacts, the MC binning can be broadened to a limit given by the theorem of Nyquist-Shannon ( $< λ / 2$ ). Although we have chosen here to model the focusing of a collimated beam with a uniform profile in the angular domain, more intricate beam profiles (Gaussian, Bessel, etc.,) can be easily implemented.²⁹

Likewise, we can conclude the following from our proof of concept in vacuo and in absorbing media: (1) the “direct extinction method” using inhomogeneous plane waves yields a very good approximation when treating problems involving absorption; (2) the assumptions adopted in the FDTD simulations (as explained in Sec. 3.2) work equally reasonably well as shown by the very good agreement with the MC results despite the numerical challenges involved. Certainly, including Fresnel reflection phenomena and extending our model to the 3-D case would be of particular importance in a next step.

Another significant improvement would be the treatment of scattering processes within the medium in which the beam is focused. A preliminary investigation has been undertaken here by calculating the direction of the energy flow in the vicinity of the focal region. This highlights a meaningful advantage of our method: We have shown that our alternative MC approach can also account for waves that show longitudinal components (nonTEM, as for example in Fig. 8). The modeling of such longitudinal components is part of the limitations of the existing vectorial MC programs: this necessitates the development of new adequate models, as initiated by Azzam.³³

Acknowledgments

This research has been supported by the Evangelisches Studienwerk e.V. Villigst ( www.evstudienwerk.de/english.html), the Baden-Württemberg Stiftung ( www.bwstiftung.de/en) and the Swiss National Science Foundation ( www.snf.ch/en), Grant No. PBBEP2_142142. We also wish to thank Joachim Wiest, Philipp Krauter, Florian Voit, Dominik Reitzle, Emanuel Simon, and Detlef Russ for fruitful discussions.

References

1.

M. Gu, “Imaging with a high numerical aperture objective,” Advanced Optical Imaging Theory, Springer, Heidelberg (2000). Google Scholar

2.

J. Ploem, “Laser scanning fluorescence microscopy,” Appl. Opt., 26 (16), 3226 –3231 (1987). http://dx.doi.org/10.1364/AO.26.003226 APOPAI 0003-6935 Google Scholar

3.

W. DenkJ. StricklerW. Webb, “Two-photon laser scanning fluorescence microscopy,” Science, 248 (4951), 73 –76 (1990). http://dx.doi.org/10.1126/science.2321027 SCIEAS 0036-8075 Google Scholar

4.

H. KanoS. MizuguchiS. Kawata, “Excitation of surface-plasmon polaritons by a focused laser beam,” J. Opt. Soc. Am. B: Opt. Phys., 15 (4), 1381 –1386 (1998). http://dx.doi.org/10.1364/JOSAB.15.001381 JOBPDE 0740-3224 Google Scholar

5.

W. ZipfelR. WilliamsW. Webb, “Nonlinear magic: multiphoton microscopy in the biosciences,” Nat. Biotechnol., 21 (11), 1369 –1377 (2003). http://dx.doi.org/10.1038/nbt899 NABIF9 1087-0156 Google Scholar

6.

N. Laiet al., “Controlling aspect ratio of focal spots of high numerical aperture objective lens in multi-photon absorption process,” Opt. Commun., 258 (2), 97 –102 (2006). http://dx.doi.org/10.1016/j.optcom.2005.07.071 OPCOB8 0030-4018 Google Scholar

7.

J. Lekner, “Polarization of tightly focused laser beams,” J. Opt. A: Pure Appl. Opt., 5 (1), 6 –14 (2003). http://dx.doi.org/10.1088/1464-4258/5/1/302 JOAOF8 1464-4258 Google Scholar

8.

H. KangB. JiaM. Gu, “Polarization characterization in the focal volume of high numerical aperture objectives,” Opt. Express, 18 (10), 10813 –10821 (2010). http://dx.doi.org/10.1364/OE.18.010813 OPEXFF 1094-4087 Google Scholar

9.

A. TafloveS. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed.Artech House, Boston (2005). Google Scholar

10.

L. WangS. L. JacquesL. Zheng, “MCml-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Prog. Biomed., 47 (2), 131 –146 (1995). http://dx.doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607 Google Scholar

11.

E. A. SergeevaA. R. KatichevM. Y. Kirillin, “Two-photon fluorescence microscopy signal formation in highly scattering media: theoretical and numerical simulation,” Quantum Electron., 40 (12), 1053 (2010). http://dx.doi.org/10.1070/QE2010v040n12ABEH014506 QUELEZ 1063-7818 Google Scholar

12.

C. K. HayakawaE. O. PotmaV. 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). http://dx.doi.org/10.1364/BOE.2.000278 BOEICL 2156-7085 Google Scholar

13.

F. CaiJ. YuS. He, “Vectorial electric field Monte Carlo simulations for focused laser beams (800 nm-2220 nm) in a biological sample,” Prog. Electromagn. Res., 142 667 –681 (2013). http://dx.doi.org/10.2528/PIER13080705 PELREX 1043-626X Google Scholar

14.

F. CaiS. 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 (1), 011022 (2014). http://dx.doi.org/10.1117/1.JBO.19.1.011022 JBOPFO 1083-3668 Google Scholar

15.

A. J. WelchM. J. C. van Gemert, Optical-Thermal Response of Laser-Irradiated Tissue, Springer, Dordrecht, Heidelberg, London, New York (2011). Google Scholar

16.

B. RichardsE. Wolf, “Electromagnetic diffraction in optical systems. ii. structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A A, 253 358 –379 (1959). http://dx.doi.org/10.1098/rspa.1959.0200 PRLAAZ 0080-4630 Google Scholar

17.

N. MetropolisS. Ulam, “The Monte Carlo method,” J. Am. Stat. Assoc., 44 (247), 335 –341 (1949). http://dx.doi.org/10.1080/01621459.1949.10483310 JSTNAL 0003-1291 Google Scholar

18.

A. LiemertA. Kienle, “Analytical solution of the radiative transfer equation for infinite-space fluence,” Phys. Rev. A, 83 (1), 015804 (2011). http://dx.doi.org/10.1103/PhysRevA.83.015804 PLRAAN 1050-2947 Google Scholar

19.

A. LiemertA. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep., 3 2018 (2013). http://dx.doi.org/10.1038/srep02018 SRCEC3 2045-2322 Google Scholar

20.

B. C. WilsonG. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys., 10 (6), 824 –830 (1983). http://dx.doi.org/10.1118/1.595361 MPHYA6 0094-2405 Google Scholar

21.

J. J. Stamnes, “Focusing of two-dimensional waves,” J. Opt. Soc. Am., 71 (1), 15 –31 (1981). http://dx.doi.org/10.1364/JOSA.71.000015 JOSAAH 0030-3941 Google Scholar

22.

J. W. Goodman, Introduction to Fourier Optics, 3rd ed.Roberts and Company Publishers, Englewood, Colorado (2005). Google Scholar

23.

M. BornE. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, Cambridge University Press, Cambridge (1999). Google Scholar

24.

A. ElmakliziJ. SchäferA. Kienle, “Simulating the scanning of a focused beam through scattering media using a numerical solution of Maxwell’s equations,” J. Biomed. Opt., 19 (7), 071404 (2014). http://dx.doi.org/10.1117/1.JBO.19.7.071404 JBOPFO 1083-3668 Google Scholar

25.

P. Debye, “Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie,” Ann. Phys., 335 (14), 755 –776 (1909). http://dx.doi.org/10.1002/(ISSN)1521-3889 ANPYA2 0003-3804 Google Scholar

26.

S. Quabiset al., “Focusing light to a tighter spot,” Opt. Commun., 179 (1–6), 1 –7 (2000). http://dx.doi.org/10.1016/S0030-4018(99)00729-4 OPCOB8 0030-4018 Google Scholar

27.

M. MatsumotoT. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul., 8 3 –30 (1998). http://dx.doi.org/10.1145/272991.272995 ATMCEZ 1049-3301 Google Scholar

28.

İ. R. ÇapoğluA. TafloveV. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express, 16 (23), 19208 –19220 (2008). http://dx.doi.org/10.1364/OE.16.019208 OPEXFF 1094-4087 Google Scholar

29.

İ. R. ÇapoğluA. TafloveV. Backman, “Computation of tightly-focused laser beams in the FDTD method,” Opt. Express, 21 (1), 87 –101 (2013). http://dx.doi.org/10.1364/OE.21.000087 OPEXFF 1094-4087 Google Scholar

30.

A. BoivinJ. DowE. Wolf, “Energy flow in the neighborhood of the focus of a coherent beam,” J. Opt. Soc. Am., 57 (10), 1171 –1175 (1967). http://dx.doi.org/10.1364/JOSA.57.001171 JOSAAH 0030-3941 Google Scholar

31.

L. CicchitelliH. HoraR. Postle, “Longitudinal field components for laser beams in vacuum,” Phys. Rev. A, 41 (7), 3727 –3732 (1990). http://dx.doi.org/10.1103/PhysRevA.41.3727 PLRAAN 1050-2947 Google Scholar

32.

R. DornS. QuabisG. Leuchs, “Sharper focus for a radially polarized light beam,” Phys. Rev. Lett., 91 (23), 233901 (2003). http://dx.doi.org/10.1103/PhysRevLett.91.233901 PRLTAO 0031-9007 Google Scholar

33.

R. M. A. Azzam, “Three-dimensional polarization states of monochromatic light fields,” J. Opt. Soc. Am. A, 28 (11), 2279 –2283 (2011). http://dx.doi.org/10.1364/JOSAA.28.002279 JOAOD6 0740-3232 Google Scholar

Biography

Arnd R. Brandes completed his bachelor’s degree in information technology in 2009 and his master’s degree in biomedical engineering in 2010. Since then, he has been working as a PhD student in the Material Optics and Imaging Group under the supervision of Alwin Kienle.

Ahmed Elmaklizi is a PhD student at the Institut für Lasertechnologien in der Medizin and Meßtechnik (ILM), Ulm, Germany, and part of the Materials Optics and Imaging Group at ILM. He received his BSc degree in electrical engineering from the German University in Cairo (GUC), and received his MS degree in communication engineering from the University of Ulm.

H. Günhan Akarçay completed his PhD degree at the Institute of Applied Physics (IAP) of the University of Bern before carrying out postdoctoral research in Montreal and Ulm. His work focuses on the modeling of light propagation in soft matter and his interests cover radiative transfer, electromagnetism, statistical physics, and biophysics. Currently, he works as a postdoctoral fellow at the IAP.

Alwin Kienle is vice director (science) of the Institut für Lasertechnologien in der Medizin and Meßtechnik (ILM), Ulm, Germany, and head of the Materials Optics and Imaging Group at ILM. In addition, he is a professor at the University of Ulm. He studied physics and received his doctoral and habilitation degrees from the University of Ulm. As postdoc, he worked with research groups in Hamilton, Canada, and in Lausanne, Switzerland.

Citation Download Citation

Arnd R. Brandes, Ahmed Elmaklizi, H. Günhan Akarcay, and Alwin Kienle "Modeling the tight focusing of beams in absorbing media with Monte Carlo simulations," Journal of Biomedical Optics 19(11), 115003 (13 November 2014). https://doi.org/10.1117/1.JBO.19.11.115003

Published: 13 November 2014

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 1 scholarly publication.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Monte Carlo methods

Finite-difference time-domain method

Optical simulations

Absorption

Computer simulations

Refractive index

Diffraction

1.

Introduction

2.

Theory