Modeling nonlinear optical microscopy in scattering media, part II. Radiation from focal volume to far-field: tutorial

Janaka C. Ranasinghesagara; Janaka C. Ranasinghesagara; Eric O. Potma; Eric O. Potma; Vasan Venugopalan; Vasan Venugopalan

doi:10.1364/JOSAA.478713

1. INTRODUCTION

Both deterministic and stochastic approaches have been taken to model the processes involved in optical microscopy of scattering media. One class of approaches considers the medium to consist of randomly distributed scatterers and applies a stochastic model to simulate electromagnetic beam propagation in the scattering media [1–4]. Another class of approaches utilizes a 3D Green’s function with a series of convolutions to calculate forward and backward scattered fields [5–7]. While such approaches may provide qualitative agreement with experimental observations, they are unable to properly quantify field characteristics for specific scattering configurations. In contrast to these approaches, the finite difference time domain (FDTD) method has been successfully used to model the propagation of tightly focused beams in scattering media for specific scattering configurations [8–10]. However, FDTD approaches require special techniques to define optical sources that possess spatially varying intensity distributions. Moreover, FDTD methods demand significant computational resources in terms of memory and run time, which scale with the total volume of the system considered [10–13].

In this two-part tutorial, we show how existing analytic electromagnetic techniques can be synthesized into a unified, comprehensive, modeling framework for nonlinear optical (NLO) microscopy in scattering media. One of the advantages of the approach presented in these tutorials is that the computational cost does not depend on the volume of the computational domain as in FDTD, but rather on the number of scatterers and number of locations where we compute the field.

This is Part II of the tutorial. In Part I [14], we formulated the mathematical representation of the incident beam, and presented several electromagnetic methods to propagate focused fields from the lens to the focal volume. The processes discussed in Part I [14] are included in the dotted gray boxes in Fig. 1. We presented three solutions to propagate focal fields in non-scattering media based on the Debye–Wolf integral (DWI), Kirchhoff’s vector integral (KVI) theorem, and Huygens–Fresnel principle (HFP). To propagate focused fields in scattering media, we considered the computationally efficient HFP solution because the DWI solution [15] is less flexible in a densely scattering medium, and the KVI solution is computationally expensive due to partial derivative calculations. We elucidated the HFP solution and included a step-by-step guide to calculating scattered fields in the focal volume and limited our scope to spherical scatterers. Since Part II uses some of the equations of Part I [14], readers are advised to read it before proceeding.

Fig. 1. Basic microscopy system (red dashed box). The processes involved in the dotted gray boxes are discussed in Part I [14] of this tutorial. The number in front of the process indicates the related sections in this tutorial.

Download Full Size | PDF

In Part II, we determine with high spatial resolution the electric field distribution proximal to the focal volume, use those results to compute the induced NLO polarization within the sample, and describe the subsequent radiation that propagates into the far field. Each element of the process is considered, beginning with sampling of the electric field near the focal volume followed by signal generation within, and radiation from, the focal volume, and concluding with far-field detection in the microscope system as shown in Fig. 1. Note that even though the excitation pulses used in NLO microscopy are broadband, for the purpose of modeling propagation effects, we maintain a monochromatic description.

The structure of Part II is as follows: in Section 2, we show how to calculate the electric field distribution near the focal volume, and in Section 3, we use these results to compute the distribution of the resulting induced nonlinear polarization within the sample. In Section 4, we treat signal radiation from the focal volume that propagates in the direction of the detector. In Section 5, we describe a method to display far-field data and continuous propagation of NLO signals in a $4f$ system. Finally, in Section 6, we provide case studies that illustrate results obtained at each stage of this process for second harmonic generation (SHG) and two-photon excitation (TPE) microscopy.

2. FOCAL SAMPLING

In Part I [14], we outlined methods to compute the electric field distribution ${\textbf{E}}({\boldsymbol \rho})$ (Fig. 1) generated by focused beam propagation in non-scattering and scattering media. The Cartesian components of ${\textbf{E}}({\boldsymbol \rho})$ are represented by ${E_x}({\boldsymbol \rho})$, ${E_y}(\boldsymbol \rho)$, and ${E_z}(\boldsymbol \rho)$, and the origin of the Cartesian coordinate system is placed at the focal point. Proper sampler placement is critical to adequately resolve the electric field distribution within the focal volume. Here, we outline two sampling configurations that we employ to record the electric field distribution within the focal volume. These field distributions, in turn, serve as an input to calculate the subsequent radiation/emission from the focal volume.

A. Cuboidal Grid Sampler

Accurate determination of the NLO signals generated in the focal region requires adequate sampling of the field within a finite volume surrounding the focal point. To accomplish this, we can define a volumetric grid sampler composed of cuboidal-shaped voxels. A volumetric grid is a 3D grid of values organized into rows, columns, and depth stacks. Each node of the grid sampler represents a field collection point as shown in Fig. 2(a). The computational cost scales linearly with the number of voxels and is independent of their size. Therefore, for a fixed volume, the computational cost scales with ${(1/\Delta R)^3}$, where $\Delta R$ is the desired voxel resolution.

Fig. 2. Intensity maps captured with a (a) cuboidal grid sampler and (b) flat sampler with a circular aperture. Size of the cuboidal volume is ${{4}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{4}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{6}}\;{\rm{\unicode{x00B5}{\rm m}}}$. Only selected slices and lines are shown for clarity (not to scale). If a solution from the KVI is considered in continuous signal propagation in Section 5.B, additional sampling locations are required at $z^\prime = z - \Delta z$ and $z^\prime = z + \Delta z$ to compute the partial derivatives.

Download Full Size | PDF

B. Planar Sampler Perpendicular to the Optical Axis

For microscopy systems with continuous signal propagation, conservation of energy demands that the total energy passing through any plane perpendicular to the optical axis remains constant. To compute the field distribution within the focal plane, we can consider a sampler plane oriented perpendicular to the optical axis and intersecting at any desired $z$ location in the vicinity of the focal point. Such a plane, possessing a fully open aperture, would capture the complete energy content of the incident beam. Practically, it is more feasible to capture the field within a finite region surrounding the focal point. In this case, we utilize a plane grid sampler in the form of a finite circular aperture [Fig. 2(b)]. If one intends to use a solution from the KVI to propagate fields in non-scattering media discussed in Section 5.B, two additional parallel plane samplers at $z^\prime = z - \Delta z$ and $z^\prime = z + \Delta z$ will be necessary to obtain partial derivatives.

3. SIGNAL GENERATION

Once we have obtained the total electrical field distribution within the focal volume, we are in a position to compute the generation of any linear or nonlinear signal of interest. In this section, we outline the approaches we utilize to calculate the generation of linear and nonlinear signals most frequently encountered in biological microscopy.

After the focal fields are captured with a cuboidal grid sampler, we can determine the signal generation efficiency at a given grid location of a voxel. For instance, for conventional fluorescence and multi-photon microscopy, we assign a fluorescence signal generation efficiency (${\le} 1$) to each grid location. For the other microscopy techniques in this section, we assign nonlinear susceptibility tensors to each voxel. The second-order nonlinear susceptibility is a tensor of rank 3, with 27 elements. The third-order nonlinear susceptibility is a tensor of rank 4, with 81 elements. More information regarding nonlinear susceptibility tensors can be found in [16–19].

A. Conventional Fluorescence Microscopy and Multi-Photon Microscopy

In single-photon and multi-photon excitation microscopy, a fluorescent molecule is excited, after which it radiates as a harmonically oscillating electric dipole. This radiation occurs incoherently relative to other radiating dipoles. We formulate the magnitude of the dipole moment $|{\textbf{p}}(\boldsymbol \rho)|$ and its excitation and emission unit vectors (${\hat{\textbf{p}}_{{\rm{ex}}}}$, ${\hat{\textbf{p}}_{{\rm{em}}}}$) separately. ${\hat{\textbf{p}}_{{\rm{ex}}}}$ and ${\hat{\textbf{p}}_{{\rm{em}}}}$ are the randomly oriented dipole moment unit vectors of the excitation and emission dipoles, respectively.

To describe the excitation of the molecular dipole, we consider four relevant cases as depicted in Fig. 3. The unit vector of the polarized excitation electric field at the molecules is given by ${\textbf{E}}({\boldsymbol \rho})/|{\textbf{E}}({\boldsymbol \rho})|$. Cases I and III consider a molecular structure with an isotropic polarizability [20,21]. In these cases, the molecular dipole is driven with the full amplitude of the excitation field, and the direction of the excited dipole ${\hat{\textbf{p}}_{{\rm{ex}}}}$ always coincides with the direction of ${\textbf{E}}({\boldsymbol \rho})$. Cases II and IV describe molecules with a transition dipole moment directed along a main axis within the molecule. We may assume a single dipole axis, and thus a specific direction ${\hat{\textbf{p}}_{{\rm{ex}}}}$. In these cases, the dipole is excited only by the component of the electric field parallel to the dipole moment [20,21].

Fig. 3. Excitation electric field unit vector ${\textbf{E}}({\boldsymbol \rho})/|{\textbf{E}}({\boldsymbol \rho})|$, dipole moment unit vector of excited dipole ${\hat{\textbf{p}}_{{\rm{ex}}}}$, and dipole moment unit vector of emission dipole ${\hat{\textbf{p}}_{{\rm{em}}}}$ in cases I–IV.

Download Full Size | PDF

In general, a randomly oriented dipole moment unit vector can be expressed as

(1)$${\hat{\textbf{p}}_{{\rm{ex/em}}}} = \left[{\begin{array}{*{20}{c}}{\cos {\phi _{\textit{dp}}}\sin {\theta _{\textit{dp}}}}\\{\sin {\phi _{\textit{dp}}}\sin {\theta _{\textit{dp}}}}\\{\cos {\theta _{\textit{dp}}}}\end{array}} \right],$$

where ${\hat{\textbf{p}}_{{\rm{ex/em}}}}$ represents either ${\hat{\textbf{p}}_{{\rm{ex}}}}$ or ${\hat{\textbf{p}}_{{\rm{em}}}}$. Angles ${\phi _{\textit{dp}}}$ and ${\theta _{\textit{dp}}}$ can be obtained by randomly sampling the surface of the unit sphere using the following expressions [22]:

(2)$$\begin{split}{{\phi _{\textit{dp}}}}&={ 2\pi \xi},\\[-3pt]{{\theta _{\textit{dp}}}}&={ \arccos ({2\xi - 1} ),}\end{split}$$

where $\xi$ is a uniformly distributed random number in the range [0,1].

In multi-photon microscopy, the magnitude of the driven dipole moment is proportional to the intensity to the power $N$, where $N$ is the order of the multi-photon process. $N = 1$ thus represents conventional fluorescence microscopy. The dipole moment $|{{\textbf{p}}^{\textbf{F}}}(\boldsymbol \rho)|$ of a molecule at point ${\boldsymbol \rho}$ can now be written as

(3)$$|{{\textbf{p}}^{\textbf{F}}}(\boldsymbol \rho)| = \left\{{\begin{array}{*{20}{l}}{\eta |{\textbf{E}}({\boldsymbol \rho} {{)|}^2}}&{\quad {\rm{Case}}\;{\rm{I}}\;{\rm{and}}\;{\rm{III}}}\\{\eta |{\textbf{E}}({\boldsymbol \rho}) \cdot {{\hat{\textbf{p}}}_{{\rm{ex}}}}{|^2}}&{\quad {\rm{Case}}\;{\rm{II}}\;{\rm{and}}\;{\rm{IV}},}\end{array}} \right.$$

where $\eta$ represents the signal generation efficiency ($\le\! 1$). Similarly, we can write the relationship for TPE ($N = 2$) and three-photon excitation ($N = 3$) microscopy as

(4)$$|{{\textbf{p}}^{{\textbf{MP}}}}(\boldsymbol \rho)| = \left \{{\begin{array}{*{20}{l}}{\eta |{\textbf{E}}({\boldsymbol \rho} {{)|}^{2N}}}&{\quad {\rm{Case}}\;{\rm{I}}\;{\rm{and}}\;{\rm{III}}}\\{\eta |{\textbf{E}}({\boldsymbol \rho}) \cdot {{\hat{\textbf{p}}}_{{\rm{ex}}}}{|^{2N}}}&{\quad {\rm{Case}}\;{\rm{II}}\;{\rm{and}}\;{\rm{IV}}{\rm{.}}}\end{array}} \right.$$

Once the molecule is excited with a probability proportional to the polarization amplitudes of Eqs. (3) and (4), the molecule may radiate through the process of fluorescence. Fluorescence emission can also be modeled with a radiating dipole model, albeit that the wavelength of emission differs from that of the excitation field. For single-photon excitation, the emission wavelength is longer than the excitation wavelength, whereas for multi-photon excitation, the emission wavelength is shorter than the excitation wavelength.

In addition, the dipole axis for the radiating molecule, ${\hat{\textbf{p}}_{{\rm{em}}}}$, need not be identical to ${\hat{\textbf{p}}_{{\rm{ex}}}}$. If the molecule is fixed on the time scale of the fluorescence lifetime $\tau$, we may assume that ${\hat{\textbf{p}}_{{\rm{em}}}}$ has a fixed orientation (cases III and IV). However, for fluorescent probes that freely diffuse through the aqueous phase, the molecule typically displays orientational motion on the time scale of $\tau$. This effectively decouples the dipole orientation for excitation and emission, and we may assume that ${\hat{\textbf{p}}_{{\rm{em}}}}$ is randomly oriented (cases I and II) and apply Eq. (1) to find its orientation [21,23] as shown in Fig. 3.

B. Second Harmonic Generation

SHG is commonly used to image biological structures that are non-centrosymmetric, including collagen, microtubules, and myosin [24]. Second harmonic signals are generated when two photons of the same frequency $\omega$ interact with a material possessing a non-vanishing second-order nonlinear susceptibility. In such a material, the photons can combine to generate a single photon at double the optical frequency of the incident photons ($2\omega$ or half the wavelength). The nonlinear polarization density $P_l^{{\textbf{SHG}}}(\boldsymbol \rho)$ generated by the electric field distribution within the focal volume for SHG can be computed using

(5)$$P_l^{{\textbf{SHG}}}(\boldsymbol \rho) = {\varepsilon _0}\sum\limits_{m,n} \;\chi _{\textit{lmn}}^{(2)}(\boldsymbol \rho)\;{E_{\omega ,m}}(\boldsymbol \rho)\;{E_{\omega ,n}}(\boldsymbol \rho),$$

where $l,m,n = x$ or $y$ or $z$. $\chi _{\textit{lmn}}^{(2)}(\boldsymbol \rho)$ is the spatial distribution of the second-order nonlinear susceptibility tensor within the sample [19]. ${E_{\omega ,m}}(\boldsymbol \rho)$ and ${E_{\omega ,n}}(\boldsymbol \rho)$ are electric field components at location ${\boldsymbol \rho}$. ${\varepsilon _0}$ is the vacuum permittivity. SHG is a coherent technique, which means that there is a phase relation between the molecular dipoles that radiate at $2\omega$ in the focal volume.

C. Third Harmonic Generation

Third harmonic generation (THG) is used to image biological interfaces formed between structures with different (non-resonant) third-order susceptibilities [25]. Strong THG signals are obtained from structures such as cell organelles, red or white blood cells, lipid droplets, adipose tissue, myelinated axons, etc [26]. In THG, three photons of the same frequency $\omega$ interact within a material with a finite third-order susceptibility and generate a single photon at $3\omega$. The nonlinear polarization density $P_l^{{\textbf{THG}}}(\boldsymbol \rho)$ generated by the electric field distribution within the focal volume for THG oscillates at $3\omega$ and can be obtained using

(6)$$P_l^{{\textbf{THG}}}(\boldsymbol \rho) = {\varepsilon _0}\sum\limits_{m,n,q} \;\chi _{\textit{lmnq}}^{(3)}(\boldsymbol \rho)\;{E_{\omega ,m}}(\boldsymbol \rho)\;{E_{\omega ,n}}(\boldsymbol \rho)\;{E_{\omega ,q}}(\boldsymbol \rho),$$

where $l,m,n,q = x$ or $y$ or $z$. $\chi _{\textit{lmnq}}^{(3)}(\boldsymbol \rho)$ is the spatial distribution of the third-order nonlinear susceptibility tensor within the sample. Like SHG, THG is a coherent technique.

D. Second-Order Sum-Frequency Generation

Second-order sum-frequency (SFG) is a coherent technique that uses an excitation field of frequency ${\omega _{{\rm{IR}}}}$ in the mid-infrared range to drive infrared-active modes in the sample. A second probe field of frequency ${\omega _{{\rm{PR}}}}$ in the visible or near-infrared range is used to generate an upconverted signal at ${\omega _{{\rm{IR}}}} + {\omega _{{\rm{PR}}}}$. SFG microscopy allows imaging with vibrational spectroscopic contrast of biological samples. SFG depends on the second-order nonlinear susceptibility, and its polarization density is given as

(7)$$P_l^{{\textbf{SFG}}}(\boldsymbol \rho) = 2{\varepsilon _0}\sum\limits_{m,n} \;\chi _{\textit{lmn}}^{(2)}(\boldsymbol \rho)\;E_m^{{\rm{IR}}}(\boldsymbol \rho)\;E_n^{{\rm{PR}}}(\boldsymbol \rho),$$

where $E_m^{{\rm{IR}}}(\boldsymbol \rho)$ are the electric field components of the mid-infrared beam, and $E_n^{{\rm{PR}}}(\boldsymbol \rho)$ are the electric field components of the probe beam.

E. Coherent Anti-Stokes Raman Scattering

Coherent anti-Stokes Raman scattering (CARS) microscopy is useful for generating images of biological samples with contrast based on Raman-active molecular vibrations [27]. Among other applications, CARS is a popular tool for visualizing lipids in cells and tissues [28,29]. Dual-color CARS uses two excitation fields, called pump (at frequency ${\omega _{\rm{P}}}$) and Stokes (at frequency ${\omega _{\rm{S}}}$) to drive a Raman-active mode at ${\omega _{\rm{P}}} - {\omega _{\rm{S}}}$. The CARS signal derives from the third-order nonlinear susceptibility of the sample and produces a polarization in the sample that oscillates at $2{\omega _{\rm{P}}} - {\omega _{\rm{S}}}$. The CARS polarization density is given as [30]

(8)$$P_l^{{\textbf{CARS}}}(\boldsymbol \rho) = {\varepsilon _0}\sum\limits_{m,n,q} \;\chi _{\textit{lmnq}}^{(3)}(\boldsymbol \rho)\;E_m^{\rm{P}}(\boldsymbol \rho)\;E_n^{\rm{P}}(\boldsymbol \rho)\;{\left[{E_q^{\rm{S}}(\boldsymbol \rho)} \right]^*},$$

where ${[{E_q^{\rm{S}}(\boldsymbol \rho)}]^*}$ is the complex conjugate of the Stokes field. To compute $P_l^{{\textbf{CARS}}}$, we first propagate the monochromatic pump and Stokes fields independently to find the focal fields. We then take the product in Eq. (8) to determine the polarization density in focus.

The polarization density expressed in Eqs. (5)–(8) can be written as ${\textbf{P}}(\boldsymbol \rho) = N\;{\textbf{p}}(\boldsymbol \rho)$, where $N$ is the molecular number density, and ${\textbf{p}}(\boldsymbol \rho)$ is nonlinearly driven molecular polarization or dipole moment.

4. SIGNAL RADIATION

For all the nonlinear signals considered in Section 3, we can apply the dipole radiation equation (DRE) [31] to propagate the generated signal toward the far field. We consider each node within our cuboidal grid to radiate as a single dipole.

We can express the radiation emitted by a single dipole located at ${\boldsymbol \rho}$ as [31–33]

(9)$$\begin{split}{\textbf{E}}({{\textbf{r}}_c})& = - \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}\left\{\vphantom{\left({\frac{1}{{{k^2}|{\textbf{R}}{|^2}}} - \frac{i}{{k|{\textbf{R}}|}}} \right)} \hat{\textbf{r}} \times (\hat{\textbf{r}} \times {\textbf{p}}(\boldsymbol \rho)) - [3\hat{\textbf{r}}\big({\hat{\textbf{r}} \cdot {\textbf{p}}(\boldsymbol \rho)} \big) - {\textbf{p}}(\boldsymbol \rho)] \right. \\&\quad\times\left. \left({\frac{1}{{{k^2}|{\textbf{R}}{|^2}}} - \frac{i}{{k|{\textbf{R}}|}}} \right)\right\} \frac{1}{{|{\textbf{R}}|}}\exp ({ik|{\textbf{R}}|} ),\end{split}$$

where ${\textbf{R}}$ is a vector from $\boldsymbol \rho$ to ${{\textbf{r}}_c}$, and $\hat{\textbf{r}}$ is its unit vector. ${\textbf{p}}(\boldsymbol \rho)$ represents the molecular polarization or dipole moment at $\boldsymbol \rho$.

We consider the use of a lens to collect the far-field signal, which we represent as a spherical reference surface with a focal length ${f_2}$. The collector location at the spherical reference surface ${{\textbf{r}}_c}$ can be expressed as (${-}{f_2}\sin {\theta _c}\cos {\phi _c}$, ${-}{f_2}\sin {\theta _c}\sin {\phi _c}$, ${f_2}\cos {\theta _c}$). The azimuthal angle ${\phi _c}$ is measured relative to the ${+}x$ axis in the counterclockwise direction, and the polar angle ${\theta _c}$ is defined with respect to the ${+}z$ axis.

We consider two categories of signal radiation: incoherent and coherent. We apply incoherent radiation for techniques such as conventional fluorescence microscopy and multi-photon fluorescence microscopy where the dipole phases are uncorrelated. For all other microscopy techniques discussed in Section 3, we apply coherent radiation.

A. Signal Radiation in Non-scattering Media

We now apply the DRE to model signal radiation from the focal region to the far field as shown in Fig. 4. For far-field detection, it can be simplified because $|{\textbf{R}}| \gg 1/k$. The result is that the signal collected by the spherical reference surface from a single dipole located at ${\boldsymbol \rho}$ can be expressed as [20,34,35]

(10)$${{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{)|_{{\rm{1DP}}}} = - \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}\left[{\hat{\textbf{r}} \times (\hat{\textbf{r}} \times {\textbf{p}})} \right]\frac{1}{{|{\textbf{R}}|}}\,\exp\! \left({ik|{\textbf{R}}|} \right),$$

where ${\textbf{p}} = {\textbf{p}}(\boldsymbol \rho) = {[{p_x},{p_y},{p_z}]^T}$. Superscript “rsc” denotes “reference surface collection” by the collection lens, and subscript ${\rm{1DP}}$ denotes one dipole. The unit vector of $|{\textbf{R}}|$, $\hat{\textbf{r}}$ is given by

(11)$${[\hat{\textbf{r}}]^T} = \left[{\begin{array}{*{20}{c}}{{{\hat r}_x}}\\{{{\hat r}_y}}\\{{{\hat r}_z}}\end{array}} \right] = \left[{\begin{array}{*{20}{c}}{({x_c} - x)/|{\textbf{R}}|}\\{({y_c} - y)/|{\textbf{R}}|}\\{({z_c} - z)/|{\textbf{R}}|}\end{array}} \right],$$

where

Fig. 4. Signal propagation from a cuboidal grid to the spherical reference surface of the far-field collector.

Download Full Size | PDF

(12)$$|{\textbf{R}}| = {\big[{{{({x_c} - x)}^2} + {{({y_c} - y)}^2} + {{({z_c} - z)}^2}} \big]^{\frac{1}{2}}}.$$

Applying vector identity ${-}\hat{\textbf{r}} \times (\hat{\textbf{r}} \times {\textbf{p}}) = {\textbf{p}} - \hat{\textbf{r}}(\hat{\textbf{r}} \cdot {\textbf{p}})$ in Eq. (10), we can write the electric field at the collector location ${{\textbf{r}}_c}$ from the $j$th dipole at $\boldsymbol \rho$ as

(13)$$\!\!\!{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c})|_j = \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}\big[{{{\textbf{p}}_j} - {{\hat{\textbf{r}}}_j}({{\hat{\textbf{r}}}_j} \cdot {{\textbf{p}}_j})} \big]\frac{1}{{|{\textbf{R}}{|_j}}}\,\exp ({ik|{\textbf{R}}{|_j}} ),\!$$

where $\hat{\textbf{r}} \cdot {\textbf{p}} = {\hat r_x}{p_x} + {\hat r_y}{p_y} + {\hat r_z}{p_z}$.

1. Coherent Radiation

In coherent radiation, electric fields are superposed to compute the final electric field at the collector. When we expand Eq. (13) for all dipoles in volume $V$, the electric field at ${{\textbf{r}}_c}$ can be expressed as

(14)$${{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}) = \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}{\iiint _V}\left[{{\textbf{p}} - \hat{\textbf{r}}(\hat{\textbf{r}} \cdot {\textbf{p}})} \right]\frac{1}{{|{\textbf{R}}|}}\,\exp ({ik|{\textbf{R}}|} ){\rm{d}}V.$$

To compute the above integral, we apply Simpson’s 1/3 rule for a 3D volume (Appendix A) as follows:

(15)$${{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}) = \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}\sum\limits_{j = 1}^{{n_\upsilon}} {\left[{\begin{array}{*{20}{c}}{{p_x} - {{\hat r}_x}(\hat{\textbf{r}} \cdot {\textbf{p}})}\\{{p_y} - {{\hat r}_y}(\hat{\textbf{r}} \cdot {\textbf{p}})}\\{{p_z} - {{\hat r}_z}(\hat{\textbf{r}} \cdot {\textbf{p}})}\end{array}} \right]_j}\frac{{{W_j}}}{{|{\textbf{R}}{|_j}}}\exp ({ik|{\textbf{R}}{|_j}} ),$$

where ${n_\upsilon}$ represents all radiating nodes, and $j$ represents the $j$th dipole at the cuboidal nodes $m,n,q$. All $\hat{\textbf{r}}$ components in the matrix are from node $m,n,q$. $|{\textbf{R}}{|_j}$ provides the distance from $({x_m},{y_n},{z_q})$ to ${{\textbf{r}}_c}$, and all unit vector components in the matrix should be computed as emanating from $m,n,q$ to ${{\textbf{r}}_c}$. The coherent intensity can then be expressed as

(16)$$I_{{\rm{coh}}}^{{\rm{rsc}}}({{\textbf{r}}_c}) \propto |{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{)|^2}.$$

2. Incoherent Radiation

For the case of incoherent radiation, we can disregard electric field interference and consider intensity contributions from individual dipoles collected at ${{\textbf{r}}_c}$. In the case of a single dipole, we can use and write the intensity at the collector location ${{\textbf{r}}_c}$ as

(17)$$I_{{\rm{incoh}}}^{{\rm{rsc}}}({{\textbf{r}}_c})|_j \propto \Big\vert{\textbf{E}}^{{\rm{rsc}}}({{\textbf{r}}_c})|_j\Big\vert^2.$$

When we consider the contribution from all dipoles in volume $V$, incoherent intensity is formulated as

(18)$$I_{{\rm{incoh}}}^{{\rm{rsc}}}({{\textbf{r}}_c}) \propto {\iiint _V}\Big\vert{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{)|_j}{\Big\vert^2}{\rm{d}}V.$$

To calculate the integral above, we apply Simpson’s 1/3 rule for a 3D volume (Appendix A). In the case of multiple random dipole orientations, the orientationally averaged incoherent intensity can be expressed as

(19)$${I_{{\rm{incoh}}}^{{\rm{rsc}}}({{\textbf{r}}_c}) \propto \left\langle {\sum\limits_\upsilon \Big\vert{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{{)|}_j}{\Big\vert^2}{W_\upsilon}} \right\rangle .}$$

Similar to Eq. (15), $\upsilon$ represents the cuboid node $m,n,q$. ${R_\upsilon}$ provides the distance from $({x_m},{y_n},{z_q})$ to ${{\textbf{r}}_c}$.

B. Radiation of Signals in Scattering Media

We next discuss the case where scatterers are present in the volume between the focal volume and the far-field collector. Note that this section employs equations in Part I and readers are advised to read Section 5 of Part I [14] before proceeding. Consider a scatterer located at $Q$ (${x_q}$, ${y_q}$, ${z_q}$) as shown in Fig. 5. We can use the DRE to determine the incident electric field on the scatterer from a dipole located at $\boldsymbol \rho$ as [8,23,34,36]

(20)$$\begin{split}\left[{\begin{array}{*{20}{c}}{E_x^{{\rm{in}}}}\\{E_y^{{\rm{in}}}}\\{E_z^{{\rm{in}}}}\end{array}} \right]& = - \frac{{{k^2}}}{{4\pi {\varepsilon _0}}}\left\{\vphantom{\left({\frac{1}{{{k^2}|{\textbf{R}}{|^2}}} - \frac{i}{{k|{\textbf{R}}|}}} \right)}\hat{\textbf{r}} \times \big({\hat{\textbf{r}} \times {\textbf{p}}(\boldsymbol \rho)} \big) - [3\hat{\textbf{r}}(\hat{\textbf{r}} \cdot {\textbf{p}}(\boldsymbol \rho)) - {\textbf{p}}(\boldsymbol \rho)]\right.\\ &\quad\times\left. \left({\frac{1}{{{k^2}|{\textbf{R}}{|^2}}} - \frac{i}{{k|{\textbf{R}}|}}} \right) \right\} \frac{1}{{|{\textbf{R}}|}}\exp \left({ik|{\textbf{R}}|} \right),\end{split}$$

where $|{\textbf{R}}|$ is the distance from $\boldsymbol \rho$ to $Q$, and the unit vector $\hat{\textbf{r}}$ is

(21)$${[\hat{\textbf{r}}]^T} = \left[{\begin{array}{*{20}{c}}{{{\hat r}_x}}\\{{{\hat r}_y}}\\{{{\hat r}_z}}\end{array}} \right] = \left[{\begin{array}{*{20}{c}}{({x_q} - x)/|{\textbf{R}}|}\\{({y_q} - y)/|{\textbf{R}}|}\\{({z_q} - z)/|{\textbf{R}}|}\end{array}} \right].$$

The unit vector $\hat{\textbf{r}}$ represents the propagation direction, ${\hat{\textbf{u}}^q}$. We can find the $\theta$ and $\phi$ angles of the propagation vector as

Fig. 5. Signal propagation from a cuboidal node to the far-field spherical reference surface in a scattering medium.

Download Full Size | PDF

(22)$$\begin{split}\theta &={ \arccos \left({{{\hat r}_z}} \right)},\\ \phi &={ \arctan \left({\frac{{{{\hat r}_y}}}{{{{\hat r}_x}}}} \right).}\end{split}$$

To apply Eq. (63) in Part I [14] for calculating the scattered fields, we first cast the electric field in terms of parallel and perpendicular components. In the far zone, the electric field lies in a plane perpendicular to $\hat{\textbf{r}}$. We can define unit vectors ${\hat{\textbf{m}}^q}$ and ${\hat{\textbf{n}}^q}$ that are perpendicular to $\hat{\textbf{r}}$ as

(23)$$\left[{\begin{array}{*{20}{c}}{{{\hat{\textbf{m}}}^q}}\\{{{\hat{\textbf{n}}}^q}}\end{array}} \right] = \left[{\begin{array}{*{20}{c}}{\cos \theta \cos \phi}&\,\,\,{\cos \theta \sin \phi}&\,\,\,{- \sin \theta}\\{- \sin \phi}&\,\,\,{\cos \phi}&\,\,\,0\end{array}} \right].$$

The orthogonal electric field components on the plane perpendicular to $\hat{\textbf{r}}$ from the $j$th dipole can be written as

(24)$${\left[{\begin{array}{*{20}{c}}{E_\parallel ^{{\rm{in}}}}\\[6pt]{E_ \bot ^{{\rm{in}}}}\end{array}} \right]_j} = \left[{\begin{array}{*{20}{c}}{{{\hat{\textbf{m}}}^q}}\\{{{\hat{\textbf{n}}}^q}}\end{array}} \right]{\left[{\begin{array}{*{20}{c}}{E_x^{{\rm{in}}}}\\[6pt]{E_y^{{\rm{in}}}}\\[6pt]{E_z^{{\rm{in}}}}\end{array}} \right]_j},$$

where ${[E_x^{{\rm{in}}},E_y^{{\rm{in}}},E_z^{{\rm{in}}}]^T}$ is given in Eq. (20).

We next assume that a field collector or a second scatterer is located at ${{\textbf{r}}_c} = ({x_c},{y_c},{z_c})$. The components of the unit vector ${\hat{\textbf{u}}^s}$ are given by

(25)$${[{\hat{\textbf{u}}^s}]^T} = \left[{\begin{array}{*{20}{c}}{({x_c} - {x_q})/|{{\textbf{R}}_s}|}\\[3pt]{({y_c} - {y_q})/|{{\textbf{R}}_s}|}\\[3pt]{({z_c} - {z_q})/|{{\textbf{R}}_s}|}\end{array}} \right],$$

where $|{{\textbf{R}}_s}|$ is the distance from $Q$ to ${{\textbf{r}}_c}$. Note that $[E_\parallel ^{\rm{in}}, E_ \bot ^{\rm{in}}]$ in Eq. (24) is analogous to $[E_\parallel ^{\rm{in}}, E_ \bot ^{\rm{in}}]$ in Eq. (63) in Part I [14]. Furthermore, the $[{\hat{\textbf{m}}^q},{\hat{\textbf{n}}^q}]$ unit vectors in Eq. (23) are analogous to the $[{\hat{\textbf{m}}^q},{\hat{\textbf{n}}^q}]$ unit vectors in Eq. (60) in Part I. Also, ${\hat{\textbf{u}}^s}$ in Eq. (25) is analogous to ${\hat{\textbf{u}}^s}$ in Eq. (62) in Part I [14].

We can follow Eqs. (73)–(75) in Part I [14] for multiple scattering and write the scattered field of a radiating dipole in the presence of multiple scatterers as follows:

(26)$${{\textbf{E}}^{{\rm{rsc,}}\,{\rm{s}}}}({{\textbf{r}}_c}{)|_j} = {\textbf{E}}_{{\rm{pri}}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c}{)|_j} + {\textbf{E}}_{{\sec}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c}{)|_j}.$$

Using the notation in Part I [14], ${\textbf{E}}_{{\rm{pri}}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c})$ and ${\textbf{E}}_{{\sec}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c})$ can be expressed as

(27)$${{\textbf{E}}_{{\rm{pri}}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c}{{)|}_j} = \sum\limits_{q = 1}^{n{\rm Scat}} \left[{\begin{array}{*{20}{c}}{{{\hat{\textbf{m}}}^s}}\\{{{\hat{\textbf{n}}}^s}}\end{array}} \right]_a^T\left[{\frac{i}{{k{d_q}}}\exp (ik{d_q})\;{\mathbb{S}_q}} \right]}\\{\left[{\begin{array}{*{20}{c}}{{\mathbb{R}_{2q\rho}}}\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}{E_\parallel ^{\rm{in}}}\\[6pt]{E_ \bot ^{\rm{in}}}\end{array}} \right]}_j}}$$

and

(28)$$\begin{split}{{\textbf{E}}_{{\sec}}^{{\rm{rsc,}}\,{\rm{s}}}({{\textbf{r}}_c}{{)|}_j} = \sum\limits_{q = 1}^{{n{\rm Scat}}} \left\{{\sum\limits_{r(\ne q) = 1}^{{n{\rm Scat}}} \left[{\begin{array}{*{20}{c}}{{{\hat{\textbf{m}}}^s}}\\{{{\hat{\textbf{n}}}^s}}\end{array}} \right]_r^T\left[{\frac{i}{{k{d_r}}}\exp (ik{d_r})\;{\mathbb{S}_r}} \right]} \right.}\\{\left. {\left[{\begin{array}{*{20}{c}}{{\mathbb{R}_{2r\rho}}}\end{array}} \right]\left({\left[{\frac{i}{{k{d_{\textit{qr}}}}}\exp (ik{d_{\textit{qr}}})\;{\mathbb{S}_q}} \right]\left[{\begin{array}{*{20}{c}}{{\mathbb{R}_{2qr}}}\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}{E_\parallel ^{\rm{in}}}\\{E_ \bot ^{\rm{in}}}\end{array}} \right]}_j}} \right)} \vphantom{\sum\limits_{q = 1}^{{n{\rm Scat}}}}\right\},}\end{split}$$

respectively. To determine ${\hat{\textbf{m}}^s}$ and ${\hat{\textbf{n}}^s}$, we use the relationships expressed in Eq. (64) in Part I [14].

The total electric field at the collector is given by the superposition of the non-scattered field [Eq. (13)] and the scattered field [Eq. (26)]:

(29)$${{\textbf{E}}^{{\rm{rsc,}}\,{\rm{ns}} + {\rm{s}}}}({{\textbf{r}}_c}{)|_j} = {{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{)|_j} + {{\textbf{E}}^{{\rm{rsc,}}\,{\rm{s}}}}({{\textbf{r}}_c}{)|_j}.$$

Rigorous computation of the scattered field for non-spherical scatterers is beyond the scope of this tutorial. The consideration of non-spherical scatterers introduces off-diagonal elements of the scattering matrix shown in Part I, Eq. (65). Moreover, the elements of the scattering matrix become functions of both polar and azimuthal angles, as opposed to the polar angle only for the case of spherical scatterers. The introduction of azimuthal angle variation also complicates the axis rotations necessary for the scattering field calculations. Nevertheless, in principle, these can all be accommodated within our framework, but will require additional computation of the scattering matrix elements. Here, we have limited our consideration to spherical scatterers and apply a formulation of Mie theory that was derived for incident plane waves to obtain the scattering matrix elements [14].

When considering dipole radiation that is incident on a stationary scatterer located in the far zone, we approximate the incident field as a plane wavelet. In the case of large spherical scatterers, the beam incident upon the scatterer may deviate significantly from a plane wave. In such instances, the generalized Lorenz–Mie theory (GLMT) [37,38] should be applied to accurately model the effects of wave curvature on the scattered field. While the treatment of GLMT is outside the scope of the current work and requies a more expensive computation of the distance-dependent elements in the scattering matrix, our framework can accommodate the use of GMLT if the scattering matrix elements are computed.

1. Coherent Radiation

In the case of coherent radiation, the total field is found from the coherent addition of the contributions from all individual dipoles in volume $V$. We may thus write

(30)$$\begin{split}{{{\textbf{E}}^{{\rm{rsc,}}\,{\rm{tot}}}}({{\textbf{r}}_c}) =}&{ {\iiint _V}{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{{)|}_j} + {{\textbf{E}}^{{\rm{rsc,}}\,{\rm{s}}}}({{\textbf{r}}_c}{{)|}_j}{\rm{d}}V}\\ = &{ {\iiint _V}{{\textbf{E}}^{{\rm{rsc,}}\,{\rm{ns}} + {\rm{s}}}}({{\textbf{r}}_c}{{)|}_j}{\rm{d}}V.}\end{split}$$

To evaluate the integral, we can apply Simpson’s 1/3 rule to evaluate the integral over the 3D grid volume (Appendix A) and write the integrand in the form of ${[{E_x}({{\textbf{r}}_c}),{E_y}({{\textbf{r}}_c}),{E_z}({{\textbf{r}}_c})]^T}$ for the non-scattering field and ${[E_x^s({{\textbf{r}}_c}),E_y^s({{\textbf{r}}_c}),E_z^s({{\textbf{r}}_c})]^T}$ for the scattered field. Now, Eq. (30) can be expressed as

Fig. 6. Comparison of the coherent intensity collected using the electric field computed in Eq. (32) and FDTD in a medium with spherical scatterers of different sizes. (i) Schematic (not to scale). Total field radiation predicted by (ii) the method in Eqs. (26)–(32) and (iii) FDTD. (iv) Intensity difference between two methods as a percentage of maximum FDTD intensity; 27 dipoles are placed within the focal volume. The central dipole is placed at the origin, and each row and each column contain three dipoles. The dipoles are spaced 50 nm apart horizontally and vertically. (a) The case of no scatterers is provided as a reference. Scatterer sizes are (b) 0.5 µm, (c) 0.75 µm, (d) 1 µm, and (e) 1.25 µm. Scatterer locations in (b)–(e) are (${-}{0.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, 1 µm, 3 µm), (${-}{1.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, 4 µm) and (2 µm, 1 µm, 5 µm). A circular $x \text{-} y$ plane collector is placed at $z = 7 \;{\rm{\unicode{x00B5}{\rm m}}}$ to collect fields. The relative refractive index of all scatterers is 1.2, and the wavelength of the incident beam is 800 nm.

Download Full Size | PDF

(31)$${{\textbf{E}}^{{\rm{rsc,}}\,{\rm{tot}}}}({{\textbf{r}}_c}) = \sum\limits_{j = 1}^{{n_\upsilon}} {\left[{\begin{array}{*{20}{c}}{{E_x}({{\textbf{r}}_c}) + E_x^s({{\textbf{r}}_c})}\\{{E_y}({{\textbf{r}}_c}) + E_y^s({{\textbf{r}}_c})}\\{{E_z}({{\textbf{r}}_c}) + E_z^s({{\textbf{r}}_c})}\end{array}} \right]_j}{W_j}$$

or

(32)$${{\textbf{E}}^{{\rm{rsc,}}\,{\rm{tot}}}}({{\textbf{r}}_c}) = \sum\limits_{j = 1}^{{n_\upsilon}} {\left[{\begin{array}{*{20}{c}}{{E_x}({{\textbf{r}}_c})}\\{{E_y}({{\textbf{r}}_c})}\\{{E_z}({{\textbf{r}}_c})}\end{array}} \right]_j}{W_j} + \sum\limits_{j = 1}^{{n_\upsilon}} {\left[{\begin{array}{*{20}{c}}{E_x^s({{\textbf{r}}_c})}\\{E_y^s({{\textbf{r}}_c})}\\{E_z^s({{\textbf{r}}_c})}\end{array}} \right]_j}{W_j},$$

where $j$ represents the $j$th radiating dipole at cuboidal node $m,n,q$. ${W_j}$ is the Simpson’s weight at node $mnq$ (Appendix A). For the case of coherent radiation, we can compute the non-scattered and scattered fields either together or separately and superpose them at the end. Finally, the intensity at the spherical reference surface following collimation can be expressed as

(33)$$I_{{\rm{coh}}}^{{\rm{rsc,}}\,{\rm{tot}}}({{\textbf{r}}_c}) \propto |{{\textbf{E}}^{{\rm{rsc,}}\,{\rm{tot}}}}({{\textbf{r}}_c}{)|^2}.$$

An example for the case of coherent radiation in a scattering medium is given in Fig. 6, where we have considered a simple scattering configuration and separately applied the method presented in Eqs. (26)–(32) and the FDTD method. The objective of this example is to compare methods for various scatterer sizes and to show the potential limitations of applying Mie theory to the simulation of dipole radiation incident on large spherical scatterers. The simulated configuration contains 27 $x$-polarized dipoles spaced 50 nm apart in both horizontal and vertical directions. The center dipole is located at the origin. Three scatterers are placed between the dipole assembly and the collector. Four different scatterer sizes are considered for this analysis.

We place an array of nodes with 50 nm spacing in the $x \text{-} y$ plane located at $z = \;7\;{\rm{\unicode{x00B5}{\rm m}}}$ to collect the coherent intensity. By selecting a location proximal to the focal plane to obtain the electric field, we can ignore the near- to far-field (NTFF) transformation [39,40], which is necessary to obtain far-field results using FDTD. In this way, we avoid the effect of extraneous factors that may influence FDTD accuracy. The differences between the two computational methods are calculated by direct subtraction of the FDTD results from the results obtained using Eq. (32) and normalized by the maximum FDTD intensity. It is evident from Fig. 6 that the prediction provided by Eqs. (26)–(32) shows good agreement with the FDTD results for smaller scatterers, but gradually deviates for larger scatterers. Note that the FDTD data cannot be assumed as the absolute reference because it may suffer from discretization errors and unavoidable reflections from the perfectly matched layer (PML) boundary condition. On the other hand, the Mie solution for light scattering that is considered here is derived for an incident plane wave. For larger scatterers located close to the source, this assumption no longer holds. These factors may account for the more pronounced differences that we observe when examining larger scatterers.

The results of Fig. 6 allow calculation of the normalized integrated intensity in the detector plane. Such a calculation using either FDTD simulations or Eqs. (26)–(32) reveals that scattering can produce increases in the integrated intensity. Note that both our results and FDTD simulations predict similar spatial variations in the coherent intensity. For larger scatterers, high-intensity spots are observed in the detector plane. When the scatterers are placed at specific locations very close to a large number of closely packed sources, they may function as focusing lenses, and thus concentrate intensity in certain regions [41].

2. Incoherent Radiation

When we consider radiation from a single dipole, the scattered and non-scattered fields are still combined coherently at the detector. But as mentioned before, the electric field interference is ignored among dipoles in the case of incoherent radiation. We can consider all dipoles in volume $V$ and express the intensity at the detector as

(34)$$\begin{split}{I_{{\rm{incoh}}}^{{\rm{rsc,}}\,{\rm{tot}}}({{\textbf{r}}_c})}&\propto{ {\iiint _V}\Big\vert\left({{{\textbf{E}}^{{\rm{rsc}}}}({{\textbf{r}}_c}{{)|}_j} + {{\textbf{E}}^{{\rm{rsc,}}\,{\rm{s}}}}({{\textbf{r}}_c}{{)|}_j}} \right){\Big\vert^2}\,{\rm{d}}V}\\ &\propto{ {\iiint _V}\Big\vert\left({{{\textbf{E}}^{{\rm{rsc,}}\,{\rm{ns}} + {\rm{s}}}}({{\textbf{r}}_c}{{)|}_j}} \right){\Big\vert^2}\,{\rm{d}}V}\\ &\propto{ {\iiint _V}I{|_j}\,{\rm{d}}V,}\end{split}$$

where $I{|_j}$ is the intensity at the detector contributed by radiation from the $j$th dipole. After applying Simpson’s 1/3 rule for integration over the 3D grid volume (Appendix A) and considering multiple dipole orientations in the simulation, the average incoherent intensity can be expressed as

(35)$$I_{{\rm{incoh}}}^{{\rm{rsc,}}\,{\rm{tot}}}({{\textbf{r}}_c}) \propto \left\langle {\sum\limits_{j = 1}^{{n_\upsilon}} I{|_j}\,{W_j}} \right\rangle .$$

5. FAR-FIELD DETECTION AND CONTINUOUS PROPAGATION

In Section 4, we used a condenser lens represented by a spherical reference surface for collecting far-field data. In some microscopy designs, the detection module incorporates a $4f$ system [42] to produce an image of the focal region onto the detector. We refer to light propagation through such an optical relay system as continuous propagation of the far-field signal. In this section, we discuss tools that can be used to describe the continuous propagation in a $4f$ system, i.e., light propagation from the spherical reference surface to the location of the detector. In addition, we also highlight methods that can be used for properly displaying the far-field data on a spherical collection surface in ways that facilitate their interpretation.

A. Displaying Far-Field Data on a Spherical Reference Surface

As mentioned before, the location ${{\textbf{r}}_c}$ at the spherical reference surface is given by (${-}{f_2}\sin {\theta _c}\cos {\phi _c}$, ${-}{f_2}\sin {\theta _c}\sin {\phi _c}$, ${f_2}\cos {\theta _c}$). We first discuss how these angle-dependent collection points can be distributed on the spherical reference surface to display the data in different formats.

When comparing the far-field collected signal with the incident field, it is convenient to choose uniformly distributed points on the spherical reference surface [43] that are rotated by 180° relative to the optical axis. The use of a uniform point distribution on the spherical surface enables comparison of the input and output fields. Displaying electric field components or intensity data at the surface in some existing approaches may require a representation of the spherical reference surface with triangular elements. Here, we discuss some alternate approaches to display the data.

We can display the directional dependence of far-field radiation by constructing a ($\phi$, $\theta$) mesh of distributed points at the spherical reference surface. To do so, we need to consider the azimuthal angle ($\phi$) in the range from 0° to 360° ($2\pi$) and the polar angle ($\theta$) in the range from 0° to ${\theta _{\rm{max}}}$, where ${\theta _{\rm{max}}} = \mathop {\sin}\nolimits^{- 1} ({\rm NA}/{n_{\rm{m}}})$. We then choose the $\phi$ mesh to contain repeated rows of gradually increasing $\phi$ values as shown in Fig. 7(a), and the $\theta$ mesh to contain repeated columns of gradually increasing $\theta$ values as shown in Fig. 7(b). Now, each pixel location of a ($\phi$, $\theta$) mesh contains a $\phi$ value from the $\phi$ mesh and a $\theta$ value from the $\theta$ mesh. These ($\phi$, $\theta$) mesh values are used in ${{\textbf{r}}_c}$ to capture the far-field data. In the case of mapping the field in the epi-propagating radiation, $f$ in ${{\textbf{r}}_c}$ should be replaced with ${-}f$. Figures 7(c) and 7(d) show the ($\theta$, $\phi$) distributed points on the forward hemispherical reference surface and the epi-hemispherical reference surface.

Fig. 7. ($\phi$, $\theta$) mesh and far-field radiation. (a) $\phi$ mesh and (b) $\theta$ mesh. ($\theta$, $\phi$) mesh maps on to a hemispherical reference surface for (c) forward collection and (d) epi-collection. The spacings between two consecutive $\phi$ values and $\theta$ values are ${\phi _1}$ and ${\theta _1}$, respectively. For full hemispherical collection, ${\theta _{\rm{max}}}$ in (c) and (d) is 90°.

Download Full Size | PDF

In another approach, the far-field information is displayed following refraction at the collection lens, as shown in Fig. 8(a). The electric field output ${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f})$ can be measured using a rectangular $x \text{-} y$ plane detector placed perpendicular to the optical axis. To plot the data without interpolation, the rectangular $x \text{-} y$ plane detector contains uniformly distributed rows and columns in the detector plane. However, when considering the point distribution at the spherical reference surface as either uniformly distributed on the spherical surface or the previously presented ($\phi$, $\theta$) assignment, the mapped data locations at the $x \text{-} y$ plane detector will not lie on a uniformly distributed grid. In this case, data interpolation is required to display 2D data, and the interpolation may obscure important details of the spatially varying fields. Alternatively, we can define uniformly distributed $x \text{-} y$ coordinates on the plane detector and map them back to a location ${{\textbf{r}}_c}$ on the spherical reference surface with a corresponding $(\phi ,\theta)$ location. Figure 8(b) provides a graphical representation of the lens as a spherical reference surface and the detector as an $x \text{-} y$ plane detector.

Fig. 8. (a) Detection of signal by a plane detector following lens refraction. (b) Graphical representation of the lens and detector. Each (${x_f}$, ${y_f}$) has a corresponding (${\phi _c}$, ${\theta _c}$). (c) Part of a $4f$ system and (d) signal propagation from the focal plane to the spherical reference surface (lens).

Download Full Size | PDF

Regardless of the distribution of ${{\textbf{r}}_c}$ locations used on the spherical reference surface, the output electric field ${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f})$ detected on the $x \text{-} y$ plane detector can be related to the electric far field on the reference sphere ${{\textbf{E}}^{{\rm{rsc}}}}({\theta _c},{\phi _c})$ using the following transformation [34,35,44]:

(36)$${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f}) = \sqrt {\frac{{{n_{\rm{m}}}}}{{{n_{{\rm{out}}}}}}} \;\frac{1}{{a({\theta _c})}}\mathbb{R}_c^{- 1} \cdot \mathbb{L}_c^{- 1} \cdot {\mathbb{R}_c} \cdot {{\textbf{E}}^{{\rm{rsc}}}}({\theta _c},{\phi _c}),$$

where ${n_{\rm{m}}}$ and ${n_{{\rm{out}}}}$ are the refractive indices of the medium before and after the lens, respectively. $\sqrt {\frac{{{n_{\rm{m}}}}}{{{n_{{\rm{out}}}}}}} \;\frac{1}{{a({\theta _c})}}$ reverses the radiometric effect provided in Eq. (7) of Part I [14]. The matrix ${\mathbb{R}_c}$ provides rotation around the $z$ axis, and ${\mathbb{L}_c}$ provides rotation around the axis perpendicular to the meridional plane. The transpose of ${\mathbb{R}_c}$ gives $\mathbb{R}_c^{- 1}$. We can write ${\mathbb{R}_c}$ and $\mathbb{L}_c^{- 1}$ as

(37)$${\mathbb{R}_c} = \left[{\begin{array}{*{20}{c}}{\cos {\phi _c}}&{\sin {\phi _c}}&0\\{- \sin {\phi _c}}&{\cos {\phi _c}}&0\\0&0&1\end{array}} \right]$$

and

(38)$$\mathbb{L}_c^{- 1} = \left[{\begin{array}{*{20}{c}}{\cos {\theta _c}}&0&{\sin {\theta _c}}\\0&1&0\\{- \sin {\theta _c}}&0&{\cos {\theta _c}}\end{array}} \right].$$

Note that ${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f})$ is rotated by 180° relative to the incident field, ${{\textbf{E}}^{{\rm{inc}}}}({x_f},{y_f})$. If we wish to rotate the detected field to match the orientation of the incident field, we can apply the following coordinate transformation:

(39)$$\left[{\begin{array}{*{20}{c}}{x_f^{{\rm{new}}}}\\{y_f^{{\rm{new}}}}\end{array}} \right] = \left[{\begin{array}{*{20}{c}}{- 1}&0\\0&{- 1}\end{array}} \right]\left[{\begin{array}{*{20}{c}}{{x_f}}\\{{y_f}}\end{array}} \right],$$

where ($x_f^{{\rm{new}}},y_f^{{\rm{new}}}$) are the new coordinates, so that ${{\textbf{E}}^{{\rm{out}}}}(x_f^{{\rm{new}}},y_f^{{\rm{new}}})$ represents the rotated electric field.

Note that the $\mathbb{R}_c^{- 1} \cdot \mathbb{L}_c^{- 1} \cdot {\mathbb{R}_c}$ transformation in Eq. (36) is not required for plotting intensities. We can write an expression for the intensity at the $x \text{-} y$ plane detector as

(40)$${I^{{\rm{out}}}}({x_f},{y_f}) = \frac{{{n_{\rm{m}}}}}{{{n_{{\rm{out}}}}}}\;{\left({\frac{1}{{a({\theta _c})}}} \right)^2}{I^{{\rm{rsc}}}}({\theta _c},{\phi _c}),$$

where ${I^{{\rm{rsc}}}}({\theta _c},{\phi _c})$ is provided by the expressions in Eq. (16), (19), (33), or (35).

B. Continuous Propagation

So far, the focal fields calculated in Part I [14] have been applied to generate linear or NLO signals from objects in the focal volume. Instead of detecting the field with a planar sampler after the collection lens, we may consider an additional ${{4}}f$ relay system that projects the field onto a detector. The continuous propagation of ${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f})$ through such a system can be described using different approaches. Here, we provide a description based on the current framework, which proceeds in two steps.

First, ${{\textbf{E}}^{{\rm{out}}}}({x_f},{y_f})$ is propagated towards a focal plane sampler containing a circular aperture [Fig. 8(c)]. Second, the field that reaches the sampler [Fig. 8(d)] is propagated towards the far-field lens. For the first step, we can use one of the three solution integrals we have presented in Section 4 of Part I [14]. For the second step, we can consider an electromagnetic propagation method such as KVI or HFP. Since a KVI-based solution requires calculation of partial derivatives, we utilize the HFP-based solution. Similar to obtaining the integral solutions in Secion 4 of Part I [14], we can follow Eqs. (26) and (27) in Part I [14] to define the electric field beyond the focal plane as

(41)$${{\textbf{E}}^{{\rm{rsf}}}}({{\textbf{r}}_f}) = \frac{{- ik}}{{2\pi}} \iint_A {\textbf{E}}({\boldsymbol \rho})\left({\frac{{1 + \cos \beta}}{2}} \right)\frac{1}{{|{\textbf{R}}|}}\exp (ik|{\textbf{R}}|)\,{\rm{d}}A,$$

where $A$ is the circular aperture, ${\textbf{r}}$ is the unit vector of ${\textbf{R}}$, and $\beta$ is the angle between $\hat{\textbf{u}}$ and ${\textbf{r}}$. The superscript rsf denotes the spherical reference surface of the lens. We can apply Simpson’s 1/3 rule (Appendix A) or another method for the evaluation of the 2D integral and write the electric far field as

Fig. 9. Nine different simulations illustrating the excitation amplitude (${{\rm{log}}_{10}}$) at the focal volume in a medium with different sizes of spherical scatterers at different locations using Eq. (77) in Part I [14]. (i) Scatterers are placed prior (“upstream”) to the focal plane; (ii) scatterers are placed beyond (“downstream”) the focal plane; (iii) scatterers are placed both prior to and beyond the focal plane. The scatterer sizes in each column are: (a) 0.5 µm, (b) 0.75 µm, and (c) 1 µm. $x \text{-} y$ plane slices are shown at $z = - {4.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{3}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{1.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, 0, 1.5 µm, 3 µm, and 4.5 µm. The scatterer locations are: (${-}{{2}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{5}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (1.5 µm, 1 µm, ${-}{{4}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (0.5 µm, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{3}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (2 µm, 1 µm, 5 µm), (${-}{1.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, 4 µm), and (${-}{0.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, 1 µm, 3 µm). Some scatterers are not visible because they are hidden behind adjacent slices.

Download Full Size | PDF

(42)$${{\textbf{E}}^{{\rm{rsf}}}}({{\textbf{r}}_f}) = \frac{{- ik}}{{2\pi}}\sum\limits_{m,n} {\textbf{E}}({{\boldsymbol \rho} _{\textit{mn}}})\left(\!{\frac{{1 + \cos {\beta _{\textit{mn}}}}}{2}}\! \right)\frac{{{W_{\textit{mn}}}}}{{{R_{\textit{mn}}}}}\exp\! \left({ik{R_{\textit{mn}}}} \right).$$

The accuracy of the far-field data depends on the size of the focal plane aperture. A bigger aperture provides better accuracy because it considers a larger area of the electric field distribution passing through the focal plane. In an ideal case using a fully open aperture at the focal plane, the far-field signal is the inverted image of the incident beam when Herschel’s condition $[{a(\theta) = 1)}]$ is satisfied.

6. CASE STUDY

In previous sections, we have presented methods for simulating field propagation and signal generation as relevant to microscopy experiments in scattering and non-scattering media. To illustrate the utility of these methods, we present a case study to examine and compare the effects of scattering on TPE and SHG microscopy.

A. Effect of Scatterers on the Excitation Field

The first step of the process is to model the excitation field. We consider an $x$-polarized incident beam with $\lambda = \;800\;{\rm{nm}}$. The focal length and NA of the lens are 1000 µm and 0.866, respectively. The refractive index of the medium is 1.333.

To analyze the focal field, we place a ${{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{9}}\;{\rm{\unicode{x00B5}{\rm m}}}$ cuboid grid detector with cubic voxels of 50 nm centered within the focal volume. The focal point is considered the origin and has a node. In this example, the scattering medium consists of six spherical scatterers at fixed locations: (${-}{{2}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{5}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (1.5 µm, 1 µm, ${-}{{4}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (0.5 µm, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{3}}\;{\rm{\unicode{x00B5}{\rm m}}}$), (2 µm, 1 µm, 5 µm), (${-}{1.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, ${-}{{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$, 4 µm), and (${-}{0.5}\;{\rm{\unicode{x00B5}{\rm m}}}$, 1 µm, 3 µm). Three scatterer diameters are considered: 0.5 µm, 0.75 µm, and 1 µm. The refractive index of the scatterers is chosen to be 1.49 providing a relative refractive index of 1.12. In this study, we consider three scenarios: (i) scatterers placed “upstream” from the focal plane, (ii) scatterers placed beyond the focal plane, and (iii) scatterers placed on both sides of the focal plane. Examination of these scenarios allows the evaluation of the effect of the scatterer location on the focal fields and the collected far-field signal.

We utilize our computational framework based on Mie scattering and apply Eq. (77) in Part I [14] to determine the excitation field. Figure 9 shows the excitation amplitude obtained for three scenarios with three different scatterer sizes. Only a few selected slices are shown in the figure for clarity. The results are normalized relative to the maximum amplitude in a non-scattering medium. To better reveal the effect of scatterers on the electric field distribution in the focal plane, in Fig. 10, we display the percent amplitude difference relative to the unscattered case and apply Eq. (B2) in Appendix B to obtain the phase difference at the focal plane relative to the non-scattering case. As shown in Fig. 10(i), for scatterers placed prior to the focal plane, the amplitude and phase distortion increases with the scatterer size. By contrast, as shown in Fig. 10(ii), scatterers placed beyond the focal field, provide minimal amplitude distortion in the focal plane. This is due to the low amplitude of the backscattered field for Mie scatterers. Note that when scatterers are small, the overall scattered field amplitude is weak, even though the portion that is backscattered is larger as compared to larger scatterers. In contrast, for large scatterers, even when the amplitude of the overall scattered field is strong, the relative amount of backward scattered light is smaller than for small scatterers. When scatterers are placed on both sides of the focal plane, the amplitude and phase distortions at the focal plane are governed primarily by the scatterers placed “upstream” from the focal plane.

Fig. 10. Percentage of amplitude difference relative to the maximum amplitude in the non-scattering case (% Amp. Diff.) and phase difference (Ph. Diff.) at the focal plane ($z = 0$) in a medium with spherical scatterers. Sizes of the scatterers are (a) 0.5 µm, (b) 0.75 µm, and (c) 1 µm. (i) Scatterers are placed prior to the focal plane; (ii) scatterers are placed beyond the focal plane; (iii) scatterers are placed on each side of the focal plane.

Download Full Size | PDF

B. Signal Generation for TPE and SHG Microscopy

With the focal field defined, we are now in a position to calculate the NLO signal. To do this, we first define an object polarized by the focal field.

For the TPE fluorescence simulation, we consider a 1 µm diameter spherical fluorescent particle, positioned at the origin of the focal volume. For simplicity, we assume that the refractive index of the fluorescence particle is the same as the surrounding medium to avoid the effect of scattering by the fluorescent particle. The excitation amplitude at the particle in a non-scattering medium is shown in Fig. 11(a). The excitation field drives molecular dipoles within the spherical fluorescence particle, which subsequently radiate incoherently. We apply Eqs. (1)–(4) to compute them for cases I and III. While the excitation wavelength is set at 800 nm, the emission wavelength is assumed to be $\lambda = 500\;{\rm{nm}}$. The 1 µm diameter spherical fluorescence particle is represented by ${{50}}\;{\rm{nm}} \times {{50}}\;{\rm{nm}} \times {{50}}\;{\rm{nm}}$ voxels with each node representing a radiating dipole.

Fig. 11. TPE and SHG excitation and far-field radiation. Left: excitation amplitude (${{\rm{log}}_{10}}$) at the focal volume for $x$-polarized incident light. Middle: far-field radiation profiles in TPE microscopy (case I with 15 random orientations of ${\hat{\textbf{p}}_{{\rm{em}}}}$ and case III) (forward only) and SHG (forward and epi). Right: far-field $x{-}y$ profile after the collection lens with ${\rm{NA}} = {0.866}$. (i) 1 µm spherical fluorescence particle; (ii) ${{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$ slab centered at the focal point (not to scale). The number above the plot represents the normalized integrated intensity.

Download Full Size | PDF

For the SHG simulation, we consider a ${{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$ slab possessing a non-vanishing second-order nonlinear susceptibility that mimics the properties of collagen fibrils aligned in the $x \text{-} y$ plane. The second-order nonlinear susceptibility tensor elements of the slab are chosen as $\chi _{\textit{xxx}}^{(2)} = 1$ and $\chi _{\textit{xyy}}^{(2)} = \chi _{\textit{xzz}}^{(2)} = \chi _{\textit{yxy}}^{(2)} = \chi _{\textit{zxz}}^{(2)} = \chi _{\textit{yyx}}^{(2)} = \chi _{\textit{zzx}}^{(2)} = 0.536$ [45]. The excitation amplitude at the slab in a non-scattering medium is shown in Fig. 11(b). We apply Eq. (5) to compute the nonlinear polarization density. Similar to TPE microscopy, the slab volume is represented by ${{50}}\;{\rm{nm}} \times {{50}}\;{\rm{nm}} \times {{50}}\;{\rm{nm}}$ voxels, and each node represents a radiating dipole. The SHG emission wavelength is 400 nm.

We use the same source objects for the simulation in a scattering medium. Similar to the non-scattering case, we apply Eqs. (1)–(4) to compute the dipole moment and the nonlinear polarization density, followed by Eq. (35) to calculate the incoherent radiation in TPE microscopy and Eq. (33) to compute the coherent radiation in SHG microscopy.

C. Far-Field TPE and SHG Signals

To visualize the far-field emission, we utilize the techniques discussed in Section 5. We place a full hemispherical ($\theta$, $\phi$) collector with ${\theta _{\rm{max}}}={ 90^ \circ}$ and radius $f$ (${{{=}1000}}\;{\rm{\unicode{x00B5}{\rm m}}}$) to visualize the entire radiation pattern in the forward hemisphere. This provides full visualization of the far-field radiation profile for both NLO modalities. We can also place a collection lens with ${\rm{NA}} = {0.866}$ and an $x \text{-} y$ plane detector to display the signal following capture and collimation with a lens. We consider Herschel’s condition ($a({\theta _c}) = 1$).

Fig. 12. Forward detected TPE signal (case I with 15 random orientations of ${\hat{\textbf{p}}_{{\rm{em}}}}$) and SHG signal in the far field. (i) Scatterers are placed prior to the focal plane; (ii) scatterers are placed beyond the focal plane; (iii) scatterers are placed both prior to and beyond the focal plane. Scatterer diameters are (a) 0.5 µm, (b) 0.75 µm, and (c) 1 µm. The number above each plot represents the integrated intensity relative to the forward non-scattering integrated intensity.

Download Full Size | PDF

Fig. 13. Epi-detected TPE signal (case I with 15 random orientations of ${\hat{\textbf{p}}_{{\rm{em}}}}$) and SHG signal in the far field. (i) Scatterers are placed prior to the focal plane; (ii) scatterers are placed beyond the focal plane; (iii) scatterers are placed both prior to and beyond the focal plane. Scatterer diameters are (a) 0.5 µm, (b) 0.75 µm, and (c) 1 µm. The number above each plot represents the integrated intensity relative to the forward non-scattering integrated intensity.

Download Full Size | PDF

After signal generation, the NLO signal propagates in a non-scattering medium and reaches the $x \text{-} y$ plane detector. Figure 11 depicts the far-field radiation profile and the $x \text{-} y$ plane detector signal. In TPE microscopy, the molecules in cases I and III are excited with the full amplitude of the excitation field because ${\hat{\textbf{p}}_{{\rm{ex}}}}$ aligns with the excitation field as shown in Fig. 3. Case I provides an isotropic radiation profile in the far field for randomized ${\hat{\textbf{p}}_{{\rm{em}}}}$. In contrast, case III provides a Hertzian dipole type radiation profile in the far field because ${\hat{\textbf{p}}_{{\rm{em}}}}$ is fixed. In SHG, the radiation profile depends not only on the polarization of the excitation field but also on the nonlinear susceptibility tensor and the excited volume of the object. As we expected, the epi-detected radiation amplitude is heavily attenuated.

Figure 12 shows the signals at the forward $x \text{-} y$ plane detector for various scattering scenarios in TPE (case I) and SHG microscopy, while Fig. 13 depicts the corresponding signals at the epi $x \text{-} y$ plane detector. The results are normalized by the maximum amplitude in a non-scattering medium. For easy comparison, the SHG results are shown under the TPE results in Figs. 12 and 13. Recall that the size of the fluorescence object used in TPE is a 1 µm sphere, and the size of the object used in SHG is a ${{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{3}}\;{\rm{\unicode{x00B5}{\rm m}}} \times {{1}}\;{\rm{\unicode{x00B5}{\rm m}}}$ slab.

In the first case (i), scatterers are placed prior to the focal plane, in which case they distort the excitation field (Figs. 9 and 10). Even though Fig. 10 shows only the percentage of amplitude change at the focal plane, it highlights the effect of scatterers on the amplitude and phase of the excitation field. In both TPE and SHG microscopy, the excitation fields are attenuated, and the degree of attenuation increases with scatterer size. In TPE, the forward radiation profile is no different from isotropic radiation because the attenuated signal incoherently radiates in a non-scattering medium. In contrast, the SHG forward signals in a non-scattering medium give rise to different spatial patterns because the generated signals are radiated coherently, and they depend on both the amplitude and phase of the excitation field.

The intensity of the maxima increases with the scatterer size. We also observe shadows of the scatterers in the far-field spatial profiles. For epi-SHG radiation, we also observe the effect of scattering on the heavily attenuated signal.

In the second case (ii), scatterers are placed beyond the focal plane. In this case, the excitation field is minimally distorted. This case bears a similarity to the epi-radiation signal in the first case. The excitation field is attenuated by scatterers, and the generated signal radiates in a scattering medium, giving rise to significant distortion of the radiated field. Note that interference-like spatial features are more prominent in the far-field SHG profiles compared to the TPE profiles. In epi-TPE radiation, the scatterers placed beyond the focal plane have no significant effect in the far field.

In the third case (iii), scatterers are placed on each side of the focal plane, in which case the results are very similar to the second case in forward radiation and to the first case in epi-radiation, because the scatterers beyond the focal plane have a dominant effect in the far field.

In conclusion, when scatterers are placed prior to the focal plane, the excitation field is distorted. As a result, SHG and TPE signals are attenuated. When scatterers are placed beyond the focal plane, the spatial profile of the generated signal is affected. When scatterers are placed both upstream and downstream relative to the focal plane, both the excitation field and the generated signal are affected. This case study illustrates that the information provided through this framework is useful for understanding the effect of scattering on microscopy signals.

D. Focus Beam Distortions

To illustrate the utility of the electromagnetic field propagation methods developed in this tutorial, we consider a sample configuration of three cells as was considered in [10].

In that work, 3D FDTD simulations were performed on cells represented by nuclei and a large number of small organelles. Each cell was represented by a cuboid with a major diameter of 15 µm and minor diameter of 13 µm. Nuclei and half of the organelles were represented by ellipsoids. The major and minor diameters of the ellipsoidal nuclei were considered as 6 µm and 5 µm, respectively. The major and minor diameters of the ellipsoidal organelles were considered as 1.5 µm and 0.5 µm, respectively. The rest of the organelles were represented by spheres with a diameter of 0.5 µm. The locations of all cellular components were chosen randomly.

We simulate a similar scenario involving three cells placed along the optical axis. We represent the nucleus as a spherical scatterer with a diameter of 5 µm and other organelles in a cell as a mixture of a total of 99 spherical scatterers with diameters of 0.5 µm, 1 µm, and 1.5 µm. Thus the three cells were represented using a total of 300 spherical scatterers. We assume that the center of each nucleus is in the $x\text{-}z$ plane and place the organelles at random locations within the cell. The refractive indices of the nuclei and organelles are 1.4 and 1.38, respectively. The medium is matched to that of cytoplasm ($n$ = 1.36 [46]), and the lens NA is 0.68 for $n=1.36$. The wavelength and filling factor [47] of the $x$-polarized Gaussian incident beam are 800 nm and 0.55, respectively. The results are shown in Fig. 14. Because the additional terms used in Gaussian source implementation and positions of the nuclei and organelles in the prior work are unknown to us, and the ellipsoids are approximated by spheres, this result does not provide an exact match with that provided in [10]. Nonetheless, there is a good qualitative agreement between the methods. Importantly, the HFP-based simulation can be completed in a small fraction of the time required to perform the FDTD simulation.

Fig. 14. Intensity (${{\rm{log}}_{10}}$) distribution in $x{-}z$ plane ($y = 0$) in the focal volume. Intensity in (a) non scattering medium and (b) a medium with three cuboidal cells that contain 300 spherical scatterers. White lines show the approximate boundaries of each cuboidal cell. Black circles represent the nuclei and organelles present at $y = 0$ slice. Other organelles in the volume are not shown.

Download Full Size | PDF

7. CONCLUSION

In Part I of the tutorial [14], we provided a framework that utilizes existing analytical electromagnetic field propagation methods to comprehensively model optical microscopy in scattering samples with fixed scatterer configurations. In Part II, we discussed methods for mapping light distributions near focus, presented signal generation within, and radiation from, the focal volume, and concluded with far-field detection in the microscope system.

Our case studies illustrate the usefulness of the information that can be provided through this framework for understanding the effect of scattering on microscopy signals. Moreover, the case studies show how this comprehensive framework can be utilized as a foundation in laser scanning microscopy. The integrated intensity of the far-field signal given in case studies represents the intensity of the signal detected with the laser placed at a fixed location within the 3D sample. By moving the source and lens laterally and performing simulations at each location, one can obtain the signal intensities required to generate a full 2D image. By setting various $z$ depths, one can generate 3D volumetric renderings.

This comprehensive framework serves as a stepping stone toward understanding factors that control the degradation of image resolution and penetration depth within scattering media in optical microscopy.

APPENDIX A: SIMPSON’S 1/3 RULE FOR NUMERICAL INTEGRATION

Simpson’s 1/3 rule is based on the use of a quadratic polynomial to approximate the function over the range of the integral.

The 1D function $g(x)$ is

(A1)$$\begin{split}\int_{{x_1}}^{{x_M}} g(x)\;{\rm{d}}x &= \frac{{{h_x}}}{3}[g({x_1}) + 4g({x_2}) + 2g({x_3}) + 4g({x_4}) \\&\quad+... + 2g({x_{M - 2}}) + 4g({x_{M - 1}}) + g({x_M}))]\\ &= \frac{{{h_x}}}{3}\sum\limits_{m = 1}^M {w_m}\;g({x_m})\\& = \sum\limits_m {W_m}\;g({x_m}),\end{split}$$

where $M$ is the number of partitions between ${x_1}$ and ${x_M}$. ${w_m}$ represents 1D Simpson’s weights. ${h_x} = ({x_M} - {x_1})/(M - 1)$, and ${W_m} = {h_x}{w_m}/3$.

The 2D function $g(x,y)$ is

(A2)$$\begin{split}\int_{{y_1}}^{{y_N}} \int_{{x_1}}^{{x_M}} g(x,y)\;{\rm{d}}x\;{\rm{d}}y& = \frac{{{h_x}}}{3}\frac{{{h_y}}}{3}\sum\limits_{m = 1}^M \sum\limits_{n = 1}^N w_m^{{\rm{2D}}}\;w_n^{{\rm{2D}}}\;g({x_m},{y_n})\\ &= \sum\limits_{m,n} {W_{\textit{mn}}}\;g({x_m},{y_n}),\end{split}$$

where $N$ is the number of partitions between ${y_1}$ and ${y_N}$. $w_m^{{\rm{2D}}}$ and $w_n^{{\rm{2D}}}$ are 2D arrays that represent 1D Simpson’s weights as shown in Fig. 15. ${h_y} = ({y_N} - {y_1})/(N - 1)$, and ${W_{\textit{mn}}} = {h_x}\;{h_y}\;w_m^{{\rm{2D}}}\;w_n^{{\rm{2D}}}/9$.

Fig. 15. $w_m^{{\rm{2D}}}$ and $w_n^{{\rm{2D}}}$ 2D arrays with 1D Simpson’s weights.

Download Full Size | PDF

The 3D function $g(x,y,z)$ is

(A3)$$\begin{split}&\int_{{z_1}}^{{z_Q}} \int_{{y_1}}^{{y_N}} \int_{{x_1}}^{{x_M}} g(x,y,z)\;{\rm{d}}x\;{\rm{d}}y\;{\rm{d}}z\\ &= \frac{{{h_x}}}{3}\frac{{{h_y}}}{3}\frac{{{h_z}}}{3}\sum\limits_{m = 1}^M \sum\limits_{n = 1}^N \sum\limits_{q = 1}^Q w_m^{{\rm{3D}}}\;w_n^{{\rm{3D}}}\;w_q^{{\rm{3D}}}\;g({x_m},{y_n},{z_q})\\ &= \sum\limits_{m,n,q} {W_{\textit{mnq}}}\;g({x_m},{y_n},{z_q}),\end{split}$$

where $Q$ is the number of partitions between ${z_1}$ and ${z_Q}$. $w_m^{{\rm{3D}}}$, $w_n^{{\rm{3D}}}$, and $w_q^{{\rm{3D}}}$ are 3D arrays that represent Simpson’s 1D weights. It is achieved by creating an $M \times N \times Q$ data cube and assigning 1D array values similar to Fig. 15. ${h_z} = ({z_Q} - {z_1})/(Q - 1)$, and ${W_{\textit{mnq}}} = {h_x}\;{h_y}\;{h_z}\;w_m^{{\rm{3D}}}\;w_n^{{\rm{3D}}}w_q^{{\rm{3D}}}/27$.

APPENDIX B: COMPLEX VECTOR PROPERTIES TO COMPUTE PHASE ANGLE

It is customary to display the phase of individual components when a complex electric field is represented by its Cartesian components. The phase of the individual components may not correctly represent the actual phase of the complex electric field. Here, we show how we apply complex vector properties [48] to compute the phase angle of a complex electric field.

Let us assume that ${\textbf{E}}$ is a complex electric field with phase angle $\alpha$, and ${\textbf{S}}$ is a reference electric field with phase angle $\beta$. We can write the following complex vector relationship [48]:

(B1)$$\begin{split}{{\textbf{E}}{{\textbf{S}}^*}}&={ |{\textbf{E}}|\exp (i\alpha)\;|{\textbf{S}}|\exp (- i\beta)}\\ &={ |{\textbf{E}}||{\textbf{S}}|\exp (i(\alpha - \beta))}\\ &={ |{\textbf{E}}||{\textbf{S}}|\exp (i\psi)}\\ &={ |{\textbf{E}}||{\textbf{S}}|(\cos \psi + i\sin \psi),}\end{split}$$

where ${{\textbf{S}}^*}$ is a conjugate of ${\textbf{S}}$. $\psi = \alpha - \beta$ gives the phase angle of ${\textbf{E}}$ relative to ${\textbf{S}}$. $\psi$ can be calculated as

(B2)$${\psi = {{\rm tan}^{- 1}}\frac{{(E_x^iS_x^r - E_x^rS_x^i) + (E_y^iS_y^r - E_y^rS_y^i) + (E_z^iS_z^r - E_z^rS_z^i)}}{{(E_x^rS_x^r + E_x^iS_x^i) + (E_y^rS_y^r + E_y^iS_y^i) + (E_z^rS_z^r + E_z^iS_z^i)}},}$$

where ${\textbf{E}} = (E_x^r + iE_x^i)\hat{\textbf{i}}+ (E_y^r + iE_y^i)\hat{\textbf{j}}+ (E_z^r + iE_z^i)\hat{\textbf{k}}$, and ${\textbf{S}} = (S_x^r + iS_x^i)\hat{\textbf{i}} + (S_y^r + iS_y^i)\hat{\textbf{j}} + (S_z^r + iS_z^i)\hat{\textbf{k}}$. Superscripts $r$ and $i$ are used to represent real and imaginary, respectively.

When ${\textbf{S}}$ is set to be a non-complex electric field, $\beta$ becomes zero, and $\psi$ provides the phase of ${\textbf{E}}$.

Funding

National Institutes of Health (R21-GM128135); National Science Foundation (CBET-1805082).

Acknowledgment

We thank Dr. Carole Hayakawa for useful discussions and some mathematical formulas of this paper. We also thank to Prof. Andrew Dunn at the University of Texas, Austin, for generously providing his FDTD source code for our use.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. X. Gan and M. Gu, “Spatial distribution of single-photon and two-photon fluorescence light in scattering media: Monte Carlo simulation,” Appl. Opt. 39, 1575–1579 (2000). [CrossRef]

2. 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, 43903 (2009). [CrossRef]

3. 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]

4. B. H. Hokr, J. N. Bixler, G. Elpers, B. Zollars, R. J. Thomas, V. V. Yakovlev, and M. O. Scully, “Modeling focusing Gaussian beams in a turbid medium with Monte Carlo simulations,” Opt. Express 23, 8699–8705 (2015). [CrossRef]

5. M. Gu, X. Fan, and C. Sheppard, “Three-dimensional coherent transfer functions in fiber-optical confocal scanning microscopes,” J. Opt. Soc. Am. A 8, 1019–1025 (1991). [CrossRef]

6. M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, “Multi-layer Born multiple-scattering model for 3D phase microscopy,” Optica 7, 394–403 (2020). [CrossRef]

7. R. Su, J. Coupland, C. Sheppard, and R. Leach, “Scattering and three-dimensional imaging in surface topography measuring interference microscopy,” J. Opt. Soc. Am. A 38, A27–A42 (2021). [CrossRef]

8. P. Török, P. R. T. Munro, and E. E. Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express 16, 507–523 (2008). [CrossRef]

9. I. R. Çapoĝlu, J. D. Rogers, A. Taflove, and V. Backman, “The microscope in a computer: image synthesis from three-dimensional full-vector solutions of Maxwell’s equations at the nanometer scale,” Prog. Opt. 57, 1–91 (2012). [CrossRef]

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

11. I. R. Çapoĝlu, A. Taflove, V. Backman, I. R. Capoglu, A. Taflove, and V. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express 16, 19208–19220 (2008). [CrossRef]

12. Z. Wu, Y. Han, J. Wang, and Z. Cui, “Generation of Bessel beam sources in FDTD,” Opt. Express 26, 28727–28737 (2018). [CrossRef]

13. D. E. McCoy, A. V. Shneidman, A. L. Davis, and J. Aizenberg, “Finite-difference time-domain (FDTD) optical simulations: a primer for the life sciences and bio-inspired engineering,” Micron 151, 103160 (2021). [CrossRef]

14. J. C. Ranasinghesagara, E. O. Potma, and V. Venugopalan, “Modeling nonlinear optical microscopy in scattering media, part I. Propagation from lens to focal volume: tutorial,” J. Opt. Soc. Am. A 40, 867–882 (2023). [CrossRef]

15. P. Török, P. D. Higdon, R. Juškaitis, and T. Wilson, “Optimising the image contrast of conventional and confocal optical microscopes imaging finite sized spherical gold scatterers,” Opt. Commun. 155, 335–341 (1998). [CrossRef]

16. R. W. Boyd, Nonlinear Optics, 3rd ed. (Academic, 2008).

17. V. V. Krishnamachari and E. O. Potma, “Focus-engineered coherent anti-Stokes Raman scattering microscopy: a numerical investigation,” J. Opt. Soc. Am. A 24, 1138–1147 (2007). [CrossRef]

18. G. de Vito, I. Tonazzini, M. Cecchini, and V. Piazza, “RP-CARS: label-free optical readout of the myelin intrinsic healthiness,” Opt. Express 22, 13733–13743 (2014). [CrossRef]

19. C. Teulon, I. Gusachenko, G. Latour, and M.-C. Schanne-Klein, “Theoretical, numerical and experimental study of geometrical parameters that affect anisotropy measurements in polarization-resolved SHG microscopy,” Opt. Express 23, 9313–9328 (2015). [CrossRef]

20. P. D. Higdon, P. Török, and T. Wilson, “Imaging properties of high aperture multiphoton fluorescence scanning optical microscopes,” J. Microsc. 193, 127–141 (1999). [CrossRef]

21. A. Diaspro, ed., Confocal and Two-Photon Microscopy: Foundations, Applications, and Advances (Wiley-Liss, 2002).

22. A. J. Welch and M. J. van Gemert, Optical-Thermal Response of Laser-Irradiated Tissue (Plenum, 1995).

23. C. Sheppard and P. Török, “An electromagnetic theory of imaging in fluorescence microscopy, and imaging in polarization fluorescence microscopy,” Bioimaging 5, 205–218 (1997). [CrossRef]

24. F. S. Pavone and P. J. Campagnola, Second Harmonic Generation Imaging (CRC Press, 2013).

25. Y. Barad, H. Eisenberg, M. Horowitz, and Y. Silberberg, “Nonlinear scanning laser microscopy by third harmonic generation,” Appl. Phys. Lett. 70, 922–924 (1997). [CrossRef]

26. D. Débarre, W. Supatto, A. M. Pena, A. Fabre, T. Tordjmann, L. Combettes, M. C. Schanne-Klein, and E. Beaurepaire, “Imaging lipid bodies in cells and tissues using third-harmonic generation microscopy,” Nat. Methods 3, 47–53 (2006). [CrossRef]

27. S. Li, Y. Li, R. Yi, L. Liu, and J. Qu, “Coherent anti-Stokes Raman scattering microscopy and its applications,” Front. Phys. 8, 1–9 (2020). [CrossRef]

28. C. L. Evans and X. S. Xie, “Coherent anti-Stokes Raman scattering microscopy: chemical imaging for biology and medicine,” Annu. Rev. Anal. Chem. 1, 883–909 (2008). [CrossRef]

29. A. Fast, J. P. Kenison, C. D. Syme, and E. O. Potma, “Surface-enhanced coherent anti-Stokes Raman imaging of lipids,” Appl. Opt. 55, 5994–6000 (2016). [CrossRef]

30. J. C. Ranasinghesagara, G. De Vito, V. Piazza, E. O. Potma, and V. Venugopalan, “Effect of scattering on coherent anti-Stokes Raman scattering (CARS) signals,” Opt. Express 25, 8638–8652 (2017). [CrossRef]

31. J. D. Jackson, Classical Electrodynamics (Wiley, 1975).

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

33. P. Török, P. R. Munro, and E. E. Kriezis, “Rigorous near-to far-field transformation for vectorial diffraction calculations and its numerical implementation,” J. Opt. Soc. Am. A 23, 713–722 (2006). [CrossRef]

34. P. Török, P. D. Higdon, and T. Wilson, “Theory for confocal and conventional microscopes imaging small dielectric scatterers,” J. Mod. Opt. 45, 1681–1698 (1998). [CrossRef]

35. M. R. Foreman and P. Török, “Computational methods in vectorial imaging,” J. Mod. Opt. 58, 339–364 (2011). [CrossRef]

36. P. R. T. Munro and P. Török, “Calculation of the image of an arbitrary vectorial electromagnetic field,” Opt. Express 15, 9293–9307 (2007). [CrossRef]

37. G. Grehan, B. Maheu, and G. Gouesbet, “Scattering of laser beams by Mie scatter centers: numerical results using a localized approximation,” Appl. Opt. 25, 3539–3548 (1986). [CrossRef]

38. J. A. Lock and G. Gouesbet, “Generalized Lorenz-Mie theory and applications,” J. Quant. Spect. Rad. Trans. 110, 800–807 (2009). [CrossRef]

39. X. Li, A. Taflove, and V. Backman, “Modified FDTD near-to-far-field transformation for improved backscattering calculation of strongly forward-scattering objects,” IEEE Antennas Wirel. Propag. Lett. 4, 35–38 (2005). [CrossRef]

40. Y. He and Q. Guo, “Accuracy improvement of near-to-far-field transformation in FDTD calculation,” in 11th European Conference on Antennas and Propagation, EUCAP (Euraap, 2017), pp. 1106–1108.

41. J. C. Ranasinghesagara, C. K. Hayakawa, M. A. Davis, A. K. Dunn, E. O. Potma, and V. Venugopalan, “Rapid computation of the amplitude and phase of tightly focused optical fields distorted by scattering particles,” J. Opt. Soc. Am. A 31, 1520–1530 (2014). [CrossRef]

42. W. Han, W. Cheng, and Q. Zhan, “Design and alignment strategies of 4f systems used in the vectorial optical field generator,” Appl. Opt. 54, 2275–2278 (2015). [CrossRef]

43. C. G. Koay, “A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere,” J. Comp. Sci. 2, 377–381 (2011). [CrossRef]

44. J. Braat and P. Török, Imaging Optics (Cambridge University, 2019).

45. C. Teulon, I. Gusachenko, G. Latour, and M.-C. Schanne-Klein, “Forward versus backward polarization-resolved SHG imaging of collagen structure in tissues,” Proc. SPIE 9712, 971218 (2016). [CrossRef]

46. A. Brunsting and P. F. Mullaney, “Differential light scattering from spherical mammalian cells,” Biophys. J. 14, 439–453 (1974). [CrossRef]

47. T. Needham, Visual Complex Analysis (Clarendon, 1999).

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

Modeling nonlinear optical microscopy in scattering media, part II. Radiation from focal volume to far-field: tutorial

Abstract

1. INTRODUCTION

2. FOCAL SAMPLING

A. Cuboidal Grid Sampler

B. Planar Sampler Perpendicular to the Optical Axis

3. SIGNAL GENERATION

A. Conventional Fluorescence Microscopy and Multi-Photon Microscopy

B. Second Harmonic Generation

C. Third Harmonic Generation

D. Second-Order Sum-Frequency Generation

E. Coherent Anti-Stokes Raman Scattering

4. SIGNAL RADIATION

A. Signal Radiation in Non-scattering Media

1. Coherent Radiation

2. Incoherent Radiation

B. Radiation of Signals in Scattering Media

1. Coherent Radiation

2. Incoherent Radiation

5. FAR-FIELD DETECTION AND CONTINUOUS PROPAGATION

A. Displaying Far-Field Data on a Spherical Reference Surface

B. Continuous Propagation

6. CASE STUDY

A. Effect of Scatterers on the Excitation Field

B. Signal Generation for TPE and SHG Microscopy

C. Far-Field TPE and SHG Signals

D. Focus Beam Distortions

7. CONCLUSION

APPENDIX A: SIMPSON’S 1/3 RULE FOR NUMERICAL INTEGRATION

APPENDIX B: COMPLEX VECTOR PROPERTIES TO COMPUTE PHASE ANGLE

Funding

Acknowledgment

Disclosures

Data availability

REFERENCES

Data availability

Cited By

Figures (15)

Equations (47)

Journal of the Optical Society of America A

Janaka C. Ranasinghesagara	https://orcid.org/0000-0002-6069-6527
Eric O. Potma	https://orcid.org/0000-0003-3916-6131
Vasan Venugopalan	https://orcid.org/0000-0003-4781-1049