Open Access
Add to BibSonomy
Abstract
The development and application of nonlinear optical (NLO) microscopy methods in biomedical research have experienced rapid growth over the past three decades. Despite the compelling power of these methods, optical scattering limits their practical use in biological tissues. This tutorial offers a model-based approach illustrating how analytical methods from classical electromagnetism can be employed to comprehensively model NLO microscopy in scattering media. In Part I, we quantitatively model focused beam propagation in non-scattering and scattering media from the lens to focal volume. In Part II, we model signal generation, radiation, and far-field detection. Moreover, we detail modeling approaches for major optical microscopy modalities including classical fluorescence, multi-photon fluorescence, second harmonic generation, and coherent anti-Stokes Raman microscopy.
© 2023 Optica Publishing Group
1. INTRODUCTION
Laser-scanning nonlinear optical (NLO) microscopy methods enable structural and functional visualization of thick cellular and tissue samples with sub-micrometer detail [1–3]. These methods rely critically on the formation of a well-defined focal spot within the tissue specimen. Light scattering in the tissue, however, distorts the amplitude and phase of the incident laser beam wavefront, which degrades the confinement of the diffraction-limited focal volume [4,5]. As a result, the focal volume within scattering samples is often expanded in both lateral and axial dimensions and generally accompanied by an overall reduction of light intensity in the focal region. This negatively impacts the imaging depth in the tissue, as well as image resolution and signal-to-noise ratio (SNR). Overcoming the effects of light scattering constitutes a grand challenge in optical microscopy with the purpose of improving imaging depth and resolution in thick tissue samples.
Counteracting the effects of scattering is challenging. Various adaptive optics approaches have been developed to compensate for wavefront distortions introduced by light scattering in the sample. A common compensation strategy is to modify the wavefront of the incident beam with a deformable mirror (DM) or spatial light modulator (SLM) to offset the optical aberrations caused by optical components and/or refractive index variations within the sample [6–8]. Optical aberrations can be largely canceled out by adjusting the incident wavefront with the conjugate of the phase aberration. Phase distortion due to the presence of discrete scatterers has been evaluated or sensed by various methods to determine the optimal phase compensation [8]. While these approaches have been successful in improving the magnitude and SNR of optical signals in laser scanning microscopy, the performance gains achieved thus far have remained rather modest [9,10]. In some cases, knowledge of the phase front is useful to compensate for scattering effects albeit at the expense of image acquisition speed. Many adaptive optics implementations for laser-scanning microscopy circumvent direct measurement of the wavefront entirely and apply phase compensation indirectly using iterative optimization of empirical parameters such as overall signal strength [11,12]. In some cases, this indirect optimization approach results in large errors [13].
The fundamental principles governing electromagnetic beam propagation and scattering are well established [14–17]. Therefore, it stands to reason that these principles can be utilized to comprehensively model and describe the effects of light scattering in laser-scanning optical microscopy. Such a model has the potential to provide fundamental insights into the factors that control the degradation of image resolution and depth within scattering media. Such a model would enable a comparative analysis of the relative sensitivity of different imaging modalities under similar scattering conditions. This approach would provide a rigorous basis for understanding the fundamental processes that contribute to wavefront distortion. This knowledge may then be leveraged to generate a priori information to dynamically synthesize wavefronts for adaptive optics applications.
Investigators have modeled various components of the optical microscopy process. One approach taken to model electromagnetic beam propagation in scattering media is to consider randomly distributed scatterers and apply a stochastic model [18–22]. Another approach is the use of convolutions with a 3D Green’s function to calculate forward and backward scattered fields [23–25]. Such approaches may provide qualitative agreement with experimental observations, but they are unable to properly account for the amplitude and phase characteristics of propagating fields for specific scatterer configurations.
Numerical Maxwell’s equations solvers such as the finite difference time domain (FDTD) or pseudo-spectral time domain (PSTD) methods are well suited for computing optical microscopy signals for specific scatterer configurations. While both FDTD and PSTD methods have been used successfully to model tightly focused beam propagation in scattering media [26–30], the process is complex and requires special techniques to define sources and minimize errors [31–34]. Moreover, these techniques require significant computational resources in terms of memory and run time, and scale with the total volume under consideration [35]. FDTD and PSTD methods also require near- to far-field (NTFF) transformation [36] to obtain the far-field results, which may be compromised by poor accuracy [37]. In contrast, the electromagnetic propagation tools that have been used successfully to model stratified media [38,39] are limited to a few specific structures possessing a high degree of symmetry [40].
Several analytical field propagation methods have been developed to calculate the effect of scattering on tightly focused laser fields [41–47]. Most methods are limited to a single scatterer placed at a specific location [41–43]. Other methods [44,45] based on the T-matrix approach [48,49] with no adaptive optics corrections are limited in scope. The method developed by Török and co-workers [46] considers multiple spherical scatterers located proximal to the focal point and applies Mie theory to compute the scattered far field. Another approach, based on Huygens–Fresnel propagators [47], shares some similarities with the scattering calculations of [46] but has the advantage of considering the effects of interparticle scattering. Using Huygens–Fresnel propagators, we have developed a computational framework that simulates focused beam propagation, followed by the generation, radiation, and detection of nonlinear microscopy signals in media containing scattering particles at specific locations [50]. We utilize elements of this framework in this tutorial.
In this tutorial, we show how existing analytical electromagnetic field propagation methods can be synthesized into an integrated, comprehensive, modeling framework for NLO microscopy in scattering media. One important advantage of the analytical modeling framework presented here is that the computational cost does not depend on the volume of the computational domain, but rather on the number of scatterers and number of locations where we compute the field. Unless one wishes to model a system with thousands of scatterers, investigators can implement and execute the methods provided in this tutorial using a desktop/laptop computer.
In Part I of this tutorial, we consider a basic optical microscope configuration with aberration-free aplanatic lenses as shown in Fig. 1. The incident beam traverses the aplanatic lens and propagates toward the focal region in either a scattering or non-scattering medium. We present the methods used to model electric field propagation. Specifically, we formulate the mathematical representation of the incident beam in Section 2 and the lens output in Section 3. In Section 4, we present and analyze several electromagnetic approaches that can be utilized to compute focal fields in non-scattering media. In Section 5, we describe methods to model the effects of discrete scatterers in the medium and compute the resultant scattered fields proximal to the focal region. In Part II of this tutorial [51], we show how to calculate the electric field distribution near the focal volume and use these results to compute the resulting induced nonlinear polarization within the sample. In addition, we treat signal radiation from the focal volume in the direction of the detector in scattering media and describe a method for displaying far-field data and continuous propagation of NLO signals in a ${4}f$ system. We conclude by providing case studies that illustrate results obtained at each stage of this modeling process for second harmonic generation (SHG) and two-photon excitation (TPE) microscopy.
Fig. 1. Excitation part of a basic microscopy system. Incident beam ${\textbf{ E}^{{\rm inc}}}({x_\infty},{y_\infty})$ is considered as a collimated beam. The number in front of the process indicates the related sections in this tutorial.
2. INCIDENT BEAM FORMULATION
We begin by considering a collimated incident beam ${\textbf{ E}^{{\rm inc}}}({x_\infty},{y_\infty})$ that propagates in a direction parallel to the optical axis. Using the Jones vector representation [52,53], the incident beam can be expressed as
A. Linearly and Circularly Polarized Beams
For an axisymmetric linearly or circularly polarized beam, we can rewrite the incident beam introduced above as [40,53]
In Eq. (2), the $\mathbb{B}\mathbb{S}$ matrix facilitates the expression of polarized illumination [54,56]:
where $c = 0$ for linearly polarized illumination and $c = i = \sqrt {- 1}$ for circularly polarized illumination. $\mathbb{P}\mathbb{V}$, defined in terms of unit intensity, can be written asB. Incorporation of a Spatial Light Modulator
We can expand this basic framework by introducing a SLM into the light path prior to the incident lens. The SLM can be used to change the phase of the transverse field before entering the imaging system. The modified incident beam can be expressed as
3. LENS AND LENS OUTPUT
We consider an aberration-free aplanatic lens that we represent geometrically using a spherical surface. We will refer to this surface as the “spherical reference surface,” and it is also known as a reference sphere [55] or a focal sphere [57]. For a flat wavefront, the radius of the spherical reference surface is ${f_1}$, which corresponds to the focal length as shown in Fig. 2. The size, curvature, and angular extent of the spherical reference surface are specified by the lens numerical aperture (NA) and focal length.
Fig. 2. Depiction of the spherical reference surface and coordinate system. The origin is placed at the focal point. We consider a refracted ray that emanates from the source and lies on a meridional plane that contains both the optical axis ($z$ axis) and refracted ray under consideration. This ray points toward the focal point and is oriented at a polar angle $\theta$ relative to the ${-}z$ axis and an azimuthal angle $\phi$ relative to the ${+}x$ axis on the plane perpendicular to the $z$ axis.
The lens output is determined by the electric field on the spherical reference surface, ${\textbf{ E}^{{\rm rs}}}(\theta ,\phi)$, which can be expressed as [53,55]
It is helpful to define various transformation matrices. We define $\mathbb{R}$ to provide azimuthal rotation around the $z$ axis [53,54]:
We define the matrix $\mathbb{L}$ to provide the rotation around the axis that lies perpendicular to the meridional plane. $\mathbb{L}$ is used to account for the rotation of the field vector, i.e., the bending of the ray toward the axis. $\mathbb{L}$ is defined as [53,54]
Equation (7) describes the effect of the lens on the incident electric field. The electric field on the spherical reference surface for a lens with a perfect anti-reflection coating (${t_\parallel} = {t_ \bot} = 1$) can be expressed as
Equation (11) provides a foundation for the focused beam propagation formulas that we provide in Section 4. Wherever ${\textbf{ E}^{{\rm rs}}}(\theta ,\phi)$ is present, it can be replaced with Eq. (11) and the corresponding ${\textbf{ E}^{{\rm inc}}}$. For linearly and circularly polarized beams, we can determine ${\textbf{ E}^{{\rm inc}}}$ using the relationships described in Eqs. (2)–(5).
4. FOCUSED BEAM PROPAGATION IN NON-SCATTERING MEDIA
In this section, we examine three methods that, based on the information provided by ${\textbf{ E}^{{\rm rs}}}(\theta ,\phi)$, can be used to compute field propagation to, and field distribution within, the focal volume for non-scattering media. The first method uses the Debye–Wolf integral (DWI), also known as the Richards and Wolf integral [40,57,62,63]. The DWI has been applied to a variety of focusing problems in homogeneous media including tightly focused fields obtained by focusing linearly polarized light [40,53–55,58]. The second approach is a solution obtained from Kirchhoff’s vector integral (KVI) [53]. The KVI theorem is a solution to the wave equation, and it is the most frequently used diffraction integral in physical optics [15,53]. The third method is based on the Huygens–Fresnel principle (HFP), which describes electric field propagation between two points [15,64]. Solutions provided by the DWI are equivalent to the Fourier spectrum of the far field [55], and solutions from the KVI and HFP are based on fundamental electromagnetic theorems or principles [15,53].
A. Theory
1. Debye–Wolf Integral
The DWI provides the electric field at any observation location ${\boldsymbol \rho}$ as follows [40,53,57]:
The infinitesimal surface area on the reference sphere, expressed in Eq. (12) in terms of the unit vectors, can be rewritten in terms of the surface angles $\theta$ and $\phi$ as follows:
By substituting Eqs. (13) and (14) into Eq. (12) and replacing $\textbf{ E}({\hat u_x},{\hat u_y})$ with ${\textbf{ E}^{{\rm rs}}}(\theta ,\phi)$, the DWI in a homogeneous medium can be formulated as [40,53,55]
Note that Eq. (15) considers the phase of ${\textbf{ E}^{{\rm rs}}}(\theta ,\phi)$ relative to the focal point. We can express the phase relative to the lens position by introducing a propagation phase $\exp (ik{f_1})$ to Eq. (15) as
Equation (16) is known as the angular representation of the focal field [55] and can be used to compute the focal electric field distribution.
2. Kirchhoff’s Vector Integral Theorem
The KVI theorem considers the electromagnetic field following its passage through an aperture. Within this formulation, this electric field is expressed as [16,53,65]
Fig. 3. Diffraction by the spherical reference surface aperture. $A$ is the spherical reference surface, and $B$ is a plane outside of the spherical reference surface whose orientation is perpendicular to the $z$ axis. $C$ is a large spherical surface centered around the observation point ${\boldsymbol \rho}(x,y,z)$. $({x_0},{y_0},{z_0})$ represents any point on the reference surface. ${ F}$ is the focal point at the origin. The origin of the Cartesian coordinate system lies at focal point ${ F}$.
After applying Kirchhoff’s boundary conditions, we can write
The partial derivative of $G$ has a negative sign because Kirchhoff’s original formulation had the normal unit vector ${\hat{\textbf{u}}}$ pointing in the ${-}z$ direction. We can safely ignore the $1/|\textbf{ r}|$ term in the bracket because the distance to the observation location from the lens $|\textbf{ r}|$ is much bigger than the wavelength ($|\textbf{ r}| \gg 1/k$, where $k = 2\pi /\lambda$). In this case, Eq. (21) can be simplified as
Substitution of Eqs. (23) and (24) into Eq. (18) gives a KVI-based solution for the propagation of the electric field toward the focal volume:
3. Huygens–Fresnel Principle
The HFP considers every point along a wavefront as a source of spherical wavelets. The electric field at any “downstream” location is given by the mutual interference of all wavelets arriving at that location. As dictated by Kirchhoff’s boundary conditions, we know surfaces $B$ and $C$ have zero contribution to the electric field propagation. We can apply the HFP for surface $A$ in Fig. 3 and write a scalar superposition integral as [64]
where ${E^{{\rm rs}}}(\theta ,\phi)$ is a scalar electric field on the spherical reference surface and $h(\textbf{ r})$ is a weighting function:The obliquity factor $K(\alpha)$ has several formulations [16,64,67–69]. The original selection of $\cos \alpha$ in Huygen’s principle [64] is an excellent choice for forward propagation because $\cos \alpha$ is non-negative for ${-}\pi /2 \le \alpha \le \pi /2$, but beyond this range, $\cos \alpha$ is negative and thus specifies backward propagating contributions that do not exist in reality for freely propagating light (Fig. 4). Using the Fresnel–Kirchhoff diffraction formula, $K(\alpha)$ can be derived as $(1 + \cos \alpha)/2$ [16].
Fig. 4. Polar plots of obliquity factors $K(\alpha) = \cos \alpha$ [64] and $K(\alpha) = (1 + \cos \alpha)/2$ [16]. While $(1 + \cos \alpha)/2$ is non-negative for the whole $\alpha$ range, $\cos \alpha$ is non-negative for only ${-}\pi /2 \le \alpha \le \pi /2$.
Even though this selection provides a non-zero contribution at $\alpha = \pi /2$, it does not provide a contribution in the backward direction. In practical focal field calculations, the choice of $\cos \alpha$ or $(1 + \cos \alpha)/2$ has no significant effect on the outcome in the limit that $\alpha$ is small. In what follows, we utilize $K(\alpha) = (1 + \cos \alpha)/2$ in Eq. (28). Using the HFP to compute the propagation of light from the lens to the focal volume, we can rewrite Eq. (28) as follows:
B. Comparison of DWI, KVI, and HFP Solution Integrals
The integral expressions obtained using these three methods have different forms. Table 1 highlights the differences in the integral expressions for the DWI, KVI, and HFP. The exponential factor in the DWI differs from the exponential factor in both the KVI and HFP. The DWI neither has an obliquity factor nor does it depend on the factor $1/|\textbf{ r}|$. Note that the DWI can describe the focal field with the same degree of accuracy as the KVI/HFP only when certain conditions are satisfied [70].
Table 1. Comparison of Solution Integrals Provided by DWI, KVI, and HFPa
A mathematical difference between the KVI and HFP is that the KVI has the partial derivative of the field. Regardless of these differences, all three methods provide the same solution at the focal point because at that location, $|\textbf{ r}| = {f_1}$, ${\psi _1} = {\psi _2}$, $\cos \alpha = 1$, and $\frac{{\partial \textbf{ E}}}{{\partial {r_0}}} = ik\textbf{ E}$. Below, we examine how the differences among these methods impact the focal field predictions at other locations.
1. Exponential Propagation Factor
The differences in the exponential propagation factors are difficult to appreciate since the form of the DWI in Eq. (16) is expressed in spherical coordinates, whereas the $|\textbf{ r}|$ term in the KVI and HFP is expressed in Cartesian coordinates. To understand the impact of the exponential factor, we compute each integral and observe the phase at the focal plane. We consider an $x$-polarized incident field, and substitute the form of $|\textbf{ r}|$ given in Eq. (22) into Eq. (25) or (29) and compare the result with the DWI expression. Figure 5(a) shows the phase of the dominant electric field component ${E_x}$ using the KVI/HFP solution and the DWI solution. In the case of the KVI/HFP solution, the phase oscillation of ${E_x}$ shows a gradual upward shift for locations away from the focal point. The greatest shift can be observed at the location farthest from the focal point. This is in contrast to the DWI prediction where the phase displays a steady oscillation without drift. The tiny spikes observed at the phase flipping locations in Fig. 5 are due to round-off and truncation errors during numerical integration.
Fig. 5. Dependency of $|\textbf{ r}|$ on the exponent. Phase of dominant component ${E_x}$ at the focal plane for $x$-polarized incidence. (a) Results obtained using the KVI/HFP solution with expression for $|\textbf{ r}|$ stated in Eq. (22) (solid blue) and DWI solution (dashed red). (b) Results obtained using KVI/HFP solution with simplified $|\textbf{ r}|$ without quadratic terms in Eq. (32) (solid blue). Incident wavelength is 0.8 µm and ${\rm NA} = {0.866}$; medium refractive index is 1.333.
This phase drift is significant within the focal plane only for locations away from the focal point. Its influence at other locations proximal to the focal point may be negligible. Nevertheless, it is instructive to analyze $|\textbf{ r}|$ in the exponent to understand its origin. Equation (22) can be rewritten as
Application of the binomial approximation and replacement of ${r_0}$ by ${f_1}$ gives
We now apply the Fraunhofer approximation, [16,53], which ignores the quadratic term in Eq. (31) when $2{f_1} \gg ({x^2} + {y^2} + {z^2})$. This provides the following approximate form for $|\textbf{ r}|$:
Utilizing spherical and cylindrical coordinates, we substitute $({x_0},{y_0},{z_0}) = ({f_1}\sin \theta \cos \phi,{f_1}\sin \theta \sin \phi,{-}{f_1}\cos \theta)$ and $(x,y,z) = (\rho \cos \varphi,\rho \sin \varphi,z)$, which allows us to write $ik|\textbf{ r}|$ in the exponent as
This argument is identical to that appearing in the exponent of the DWI in Eq. (16). The introduction of this simplified form of $|\textbf{ r}|$ in the exponent in Eqs. (25) or (29) diminishes the phase drift for locations away from the focal point on the focal plane as shown in Fig. 5(a), indicating that the quadratic terms in $\textbf{ r}$ are those that give rise to the phase drift away from the focal spot. Nevertheless, we will utilize the full expression for $|\textbf{ r}|$ as shown in the exponent of Eqs. (25) and (29) in the remainder of this tutorial.
2. Evaluation of the Partial Derivative in the KVI
The KVI contains the partial derivative of the field in Eq. (25). To compute this derivative, we need to determine the electric fields of two spherical reference surfaces with radii equal to ${f_1} + \Delta {r_0}$ and ${f_1} - \Delta {r_0}$. However, this computation can be simplified for certain wavefronts. Let us consider a monochromatic wavefront in the form of
When such a wavefront is incident upon the lens, we can write the partial derivative in Eq. (25) as a finite difference:
We see that under these conditions, the solution of the KVI is identical to the integral expression obtained within the HFP [Eq. (29)].
3. Multiplication Factors
The KVI and the HFP both contain an obliquity factor in the integrand, which is absent in the DWI. Figure 4 shows the obliquity factor variation with $\alpha$. When $\alpha$ is very small and the lens has a long focal length, the effect of the obliquity factor is negligible. But if we consider a location far away from the focal point, $\alpha$ becomes larger and the results from the KVI/HFP will deviate from the DWI. Another difference we observe in Table 1 is the distance-dependent attenuation factor. The DWI uses $1/f$, while the KVI/HFP have a $1/|\textbf{ r}|$ term. To observe the combined effect of $\alpha$ and the attenuation factor, in Fig. 6, we compute the focal field amplitude along the $x$ axis, $y$ axis, and $z$ axis using the solution from the KVI/HFP and the DWI for an incident plane wave. We also include the solution from the KVI/HFP with simplified $|\textbf{ r}|$ [Eq. (32)] in the exponent to avoid the impact from different exponential factors. The results show negligible differences between the different approaches.
Fig. 6. Computed focal field amplitudes (${{\rm log}_{10}}$) along the (a) $x$ axis, (b) $y$ axis, and (c) $z$ axis, using the solution from the DWI (dashed blue), KVI/HFP (solid red), and KVI/HFP with simplified $|\textbf{ r}|$ in the exponent (dashed red). Since the curves overlap without visible differences, only the KVI/HFP solution is fully visible. Inset on the right magnifies the dashed box located at 10 µm from the focal origin. Incident wavelength is 0.8 µm and ${\rm NA} = {0.866}$; medium refractive index is 1.333.
C. Numerical Evaluations
We have formulated the integral expressions using the DWI [Eq. (16)], KVI [Eq. (25)], and HFP [Eq. (29)] to represent the focal field. Here, we discuss how those integrals can be computed numerically.
1. Debye–Wolf Integral
Solutions to the DWI [Eq. (16)] have been obtained for different illumination conditions after substitution with the following identity [53,55,57]:
This substitution reduces the 2D integral to a 1D integral. For commonly used polarization states of the incident field, we can write the following 1D integrals [53,55,57]:
We can solve these integrals numerically and express $\textbf{ E}({\boldsymbol \rho})$ for different polarizations as linearly $x$ polarized:
2. KVI and HFP Integrals
To solve the area integral expressions for the KVI and HFP, we may express ${\rm d}A$ in either Cartesian (${\rm d}x {\rm d}y$) or spherical coordinates ($\sin \theta {\rm d}\phi {\rm d}\theta$), before proceeding with the complicated double integration. Another approach that can be applied for this integral is to use a set of uniformly distributed radiating points on the spherical reference surface, so that surface element ${\rm d}A$ can be approximated as
We use Eq. (48) to express the integral in Eqs. (25) and (29) as a summation. The solution obtained from the KVI [Eq. (25)] can now be written as
The solution obtained from the HFP [Eq. (29)] can be written as
Fig. 7. Computed electric fields in the focal volume in a non-scattering medium for $x$-polarized incident plane wave. The solution of the HFP is applied to compute (a) amplitude (${{\rm log}_{10}}$) of selected slices in the focal volume. $x - y$ plane slices are shown at $z = - {4.5}\;{\unicode{x00B5}{\rm m}}$, ${-}{3}\;{\unicode{x00B5}{\rm m}}$, ${-}{1.5}\;{\unicode{x00B5}{\rm m}}$, 0, 1.5 µm, 3 µm, and 4.5 µm. (b) Amplitude and (c) phase of electric field components at the focal plane ($z = {0}$). Incident wavelength is 0.8 µm and ${\rm NA} = {0.866}$; medium refractive index is 1.333.
D. Selection of Focused Beam Propagation Method
We have presented three integral methods that can be employed to model focused beam propagation in non-scattering media. We have provided a comparative analysis of these methods and shown how to compute them numerically. The DWI provides an efficient method for analyzing focal fields in homogeneous and stratified media [39,53,73,74] and is considered a good approximation to a well-posed boundary value problem [70]. However, any attempt to use the DWI in a scattering medium requires using either a representation of the scattering medium in $k$-space or an approach similar to that provided by Török and co-workers [46]. Representation of the scattering medium in $k$-space requires an evaluation of the effects of scattering and the propagation of the resulting $k$ vector distribution before the application of the Fourier transform. To our knowledge, such a process has not been developed and applied for computing propagation in deterministic scattering media. The approach provided by Török’s group [46] does not account for interparticle scattering and is not easily applied to a medium with multiple scatterers. By contrast, solutions provided by the KVI and HFP are in a form that is readily applicable for computing both direct scattering and inter-particle scattering necessary for modeling various microscopy modalities. The KVI solution is computationally expensive because it requires the computation of electric field gradients at the lens. However, we have shown in Eqs. (34)–(36) that for converging, diverging, and plane wavefronts, this computation can be simplified, and that the simplified version is identical to the HFP. The HFP [Eq. (29)] is computationally efficient because it does not require the computation of electric field gradients. Our prior use of the HFP method has provided results that agree well with FDTD results [47].
5. FOCUSED BEAM PROPAGATION IN SCATTERING MEDIA
In this section, we outline an approach for describing the interaction of the focused beam wavefront with scatterers found in the medium between the focusing lens and the focal volume. While many other approaches have been discussed [41–46], each method has certain limitations. Here, we describe an approach grounded in the fundamentals of scattering theory [14,15,17,47,48,75–77] and provide a comprehensive framework that can be applied to a variety of microscopy modalities under different illumination conditions.
A. Medium with a Single Scatterer
First, we discuss the modeling of focused beam propagation in the presence of a single scatterer within the medium. Second, we extend the approach for use in the case of multiple scatterers. Figure 8 depicts the graphical representation for the case of a single scatterer.
Fig. 8. Scattering by a single scatterer. (a) Complete process including the electric field at the spherical reference surface ($L$), incident electric field at the scatterer $(Q)$, and scattered field at the observation point ($O$). (b) Graphical representation of electric field components $({E_\parallel ^l,E_ \bot ^l})$ at the spherical reference surface before transformation and $({E_\parallel ^{\textit{lq}},E_ \bot ^{\textit{lq}}})$ after transformation. $E_\parallel ^l$ is parallel to the meridional plane, and $E_\parallel ^{\textit{lq}}$ is parallel to the incident plane. (c) Unit vectors of (b). (d) Graphical representation of incident field at the scatterer relative to the incident plane ($E_\parallel ^q$, $E_ \bot ^q$) and relative to the scattering plane $({E_\parallel ^{\rm{in}},E_ \bot ^{\rm{in}}})$. (e) Unit vectors of (d). (f) Graphical representation of scattered field at the observation location $O$. $E_\parallel ^{{\rm in}}$ and $E_\parallel ^s$ are parallel, and $E_ \bot ^{{\rm in}}$ and $E_ \bot ^s$ are normal to the scattering plane. (g) Unit vectors of (f). The spherical reference surface is drawn horizontally such that the scattering diagram aligns with the illustrations in other references [15,17,48]. Not to scale.
1. Electric Field Components and Unit Vectors
Up to this point, we have represented the electric fields using their Cartesian components so that the total electric field can be determined through the vector addition of each Cartesian component. However, when performing scattering calculations, we must compute the incident electric field components that lie parallel and perpendicular to a scattering plane [15,17]. As shown in Fig. 8(a), the scattering plane contains the source, a scatterer, and an observation location. We can start with the parallel and perpendicular field components at the lens and apply the required coordinate transformations to obtain the incident field at the scatterer. It is also convenient to define unit vectors ${\hat{\textbf{m}}}$ and ${\hat{\textbf{n}}}$ of a local orthonormal coordinate system to track the parallel and perpendicular electric field components as well as unit vector ${\hat{\textbf{u}}}$ that represents the propagation direction [47,75].
2. Parallel and Perpendicular Components at the Lens
Consider a collimated beam incident upon the lens at location $L = ({f_1}\sin {\theta _l}\cos {\phi _l},{f_1}\sin {\theta _l}\sin {\phi _l}, - {f_1}\cos {\theta _l})$. The first step is to consider the meridional plane as shown in Fig. 2 and compute the electric field components parallel and perpendicular to it, as follows:
To track the parallel and perpendicular components, we can determine unit vectors ${{\hat{\textbf{m}}}^l} = (\hat m_x^l,\hat m_y^l,\hat m_z^l)$ and ${{\hat{\textbf{n}}}^l} = (\hat n_x^l,\hat n_y^l,\hat n_z^l)$ at location $L$ for components $E_\parallel ^l(\theta ,{\phi _l})$ and $E_ \bot ^l({\theta _l},{\phi _l})$ in the following way:
The unit vector that represents the propagation direction ${{\hat{\textbf{u}}}^l}$ points towards the focal point $F$. It is given by (${-}\sin {\theta _l}\cos {\phi _l}$, ${-}\sin {\theta _l}\sin {\phi _l}$, $\cos {\theta _l}$).
3. Electric Field Propagation from Lens to Scatterer
Consider the placement of a scatterer at location $Q({\rho _q}\cos {\varphi _q},{\rho _q}\sin {\varphi _q},\;{z_q})$ as shown in Fig. 8(a). We aim to compute the incident electric field components $E_\parallel ^{\rm{in}}$ and $E_ \bot ^{\rm{in}}$ at location $Q$. We start with $E_\parallel ^l$ and $E_ \bot ^l$ and transform them to $E_\parallel ^{\textit{lq}}$ and $E_ \bot ^{\textit{lq}}$ at location $L$. We then propagate the electric field towards location $Q$ and apply another transformation to obtain $E_\parallel ^{\rm{in}}$ and $E_ \bot ^{\rm{in}}$.
First, we find the direction of propagation of the wavelet ${{\hat{\textbf{u}}}^{\textit{lq}}}$ that points toward $Q$ at location $L$:
We can now write the relationship between ${[E_\parallel ^{\textit{lq}},E_ \bot ^{\textit{lq}}]^T}$ and ${[E_\parallel ^l,E_ \bot ^l]^T}$ as
The unit vectors of the electric field components in Eq. (56) are shown in Fig. 8(c), and are given by
To account for the finite volume of the scatterer, we must consider separately each wavelet that reaches the scatterer. Let us consider the radiation from a single radiating point of the spherical reference surface that represents a partial field incident upon the scatterer at $Q$. We apply the solution provided by the HFP [Eq. (51)] to propagate a polarized wavelet from $L$ towards $Q$ and write the partial field ($E_\parallel ^q$, $E_ \bot ^q$) in Fig. 8(d) as
The local unit vectors at $Q$ remain unchanged. ${[{{\hat{\textbf{m}}}^q},{{\hat{\textbf{n}}}^q},{{\hat{\textbf{u}}}^q}]^T}$ is given by
We now know the partial field incident upon the scatterer at $Q$. However, its parallel and perpendicular components are given relative to the incident plane. Thus, the next step is to compute the parallel and perpendicular components relative to the scattering plane, as shown in Fig. 8(f). Similar to Eq. (55), we can compute the transformation angle pair (${\phi _s}$ and ${\theta _s}$) for this transformation as
We next transform the parallel and perpendicular components of ${E^q}$, which are defined relative to the incident plane, in terms of the parallel and perpendicular field components relative to the scattering plane, i.e., ${[E_\parallel ^{\rm{in}},E_ \bot ^{\rm{in}}]^T}$, as shown in Fig. 8(f). This transformation is given as
The subsequent step is to obtain the unit vectors of $E_\parallel ^{\rm{in}}$ and $E_ \bot ^{\rm{in}}$ as shown in Fig. 8(e). ${{\hat{\textbf{m}}}^{\rm{in}}}$ and ${{\hat{\textbf{m}}}^s}$ lie in the same plane, and ${{\hat{\textbf{n}}}^{\rm{in}}}$ is collinear with ${{\hat{\textbf{n}}}^s}$ [Fig. 8(g)]. When $|\cos {\theta _s}| = 1$, we can apply the azimuthal angle transformation given in Eq. (63) to compute unit vectors ${{\hat{\textbf{m}}}^s}={{\hat{\textbf{m}}}^{\rm{in}}} = (\hat m_x^q\cos {\phi _s},\hat m_y^q\sin {\phi _s},\hat m_z^q)$ and ${{\hat{\textbf{n}}}^s}= {{\hat{\textbf{n}}}^{\rm{in}}} = ({-}\hat n_x^q\sin {\phi _s},\hat n_y^q\cos {\phi _s},\hat n_z^q)$. When $|\cos {\theta _s}| \ne 1$, we can write orthonormal unit vector relationships to determine ${{\hat{\textbf{n}}}^s}$ and ${{\hat{\textbf{n}}}^{\rm{in}}}$ and then use that information to obtain ${{\hat{\textbf{m}}}^s}$ and ${{\hat{\textbf{m}}}^{\rm{in}}}$ [17]. This yields the following relations:
4. Computation of the Scattered Field
The rigorous computation of the scattered field is a tedious process and requires the evaluation of complex integral expressions [14,15,17] that follow from the solution for optical scattering produced by a dielectric particle. Once the incident electric field components are defined, we can consider the following simplified relationship to calculate the scattered far field [14,15,17,76]:
To represent a matrix multiplication, we can place the common multiplier inside the scattering amplitude matrix and express the partial scattered field detected at the observation location for the $j$th wavelet as
To compute the full scattered field, we need to consider wavelets emanating from all locations on the spherical reference surface. However, the partial electric fields provided by Eq. (66) cannot be directly summed because the reference frames of parallel and perpendicular components are not identical across all the partial fields considered. Superposition can be performed only once all local parallel and perpendicular components are converted to global Cartesian components. This conversion is accomplished using the following relationship:
To obtain the complete scattered field, we apply Eq. (59) for each wavelet that is oriented towards $Q$ from the spherical reference surface and use it in Eqs. (66) and (67). The total scattered field at the observation location ${\boldsymbol \rho}$ due to a single scatterer can then be expressed as
Note that Eq. (68) was derived for an observation location in the “far zone.” Here, the far zone represents locations that satisfy $d \gg {\rm max}(\lambda , a, k{a^2}/2$), where $\lambda$ is wavelength, and $a$ is the radius of the scatterer [17,76]. The distance between any point within the scatterer and the observation location should be greater than $\lambda$. For scatterers much smaller than the wavelength, the inequality $d \gg a$ defines the far-zone condition. For scatterers much larger than the wavelength, the inequality $d \gg k{a^2}/2$, which defines the Fraunhofer distance, determines the distance beyond which the far-zone expressions are valid [76].
5. Scattering Amplitude Matrix
We now take a closer look at the scattering amplitude matrix defined in Eq. (65). This matrix can be used for various scatterer shapes [48]. Here, we discuss how to obtain the scattering amplitude matrix components for spherical scatterers, which is a good approximation for the scattering objects relevant for optical microscopy of biological samples [50,80–82]. Rigorous computation of the scattered field for non-spherical scatterers is more complicated and beyond the scope of this tutorial.
For spherical scatterers, ${S_3}({\theta _s},{\phi _s}) = {S_4}({\theta _s},{\phi _s}) = 0$, and the azimuthal dependence of remaining ${S_1}({\theta _s},{\phi _s})$, ${S_2}({\theta _s},{\phi _s})$ can be ignored. In this case, we may write the scattering amplitude matrix for spherical scatterers as
6. Application of Near-Field Lorenz–Mie Theory
The scattering amplitude matrix components ${S_1}({\theta _s})$ and ${S_2}({\theta _s})$ in the Mie solution are known as asymptotic scattering amplitudes at infinity. However, when the observation location or another spherical scatterer is located at a distance that does not satisfy the far-zone criteria specified above, we can apply results from near-field Lorenz–Mie theory [15,77,83] to obtain the scattered field in lieu of using ${S_1}({\theta _s})$ and ${S_2}({\theta _s})$. Near-field Lorenz–Mie theory contains distance-dependent components that are not fully transverse. The scattered field contains a radial component (${E_r}$) in addition to polar (${E_\theta}$) and azimuthal (${E_\phi}$) components. However, we can ignore the effect of the radial component because it decays rapidly within an optical wavelength, and its effect on the overall result is negligible for well-separated scatterers. Assigning ${E_\theta}$ as the parallel component and ${E_\phi}$ as the perpendicular component, we can write a relationship similar to Eq. (65) as
Here, $\xi (kd)$ is the Riccati–Bessel function, and ${\xi ^\prime}(kd)$ is its derivative [14,15]. $\mathbb{E}(kd)$ represents a distance-dependent scattering matrix.
Note that Eq. (65) includes the phase term $\exp (ikd)$, which is not explicitly found in Eq. (71). The scattering amplitudes ${S_1}({\theta _s},{\phi _s})$, ${S_2}({\theta _s},{\phi _s})$, ${S_3}({\theta _s},{\phi _s})$, and ${S_4}({\theta _s},{\phi _s})$ in Eq. (65) are not distance-dependent functions and require a separate propagation phase to compute the scattered field at the observation location. Functions $\xi (kd)$ and ${\xi ^\prime}(kd)$ in Eq. (71), on the other hand, already include the propagation phase.
B. Medium with Multiple Scatterers
We now discuss how to extend the equations developed in Section 5.A to describe wave propagation in a medium containing multiple scatterers. When considering each scatterer independently, we can determine the scattered field due to primary scattering at the observation location. Additionally, we must consider the potential for multiple scattering, i.e., the possibility that the scattered field from one scatterer may act as an incident field on other scatterers and vice versa [Fig. 9(a)]. Fortunately, in many scenarios of practical interest, the scattered field undergoes considerable geometric attenuation before it is incident on another scatterer. Such an attenuated incident field produces an even weaker scattered field. Thus, the influence of this secondary scattered field on other scatterers, as shown in Fig. 9, is often negligibly small in many relevant cases.
Fig. 9. (a) Primary, secondary, and tertiary scattering. The scattered field of the $Q$th scatterer becomes an incident on the $R$th scatterer and creates secondary scattering. The secondary scattering reaches the observation location. In addition, it becomes an incident on the $Q$th scatterer and produces tertiary scattering. ${d_q}$ is distance from the center of the $Q$th scatterer to $\rho$, ${d_r}$ is distance from the center of the $R$th scatterer to $\rho$, and ${d_{\textit{qr}}}$ is distance from the center of the $Q$th scatterer to the center of the $R$th scatterer. (b) ${{\rm Log}_{10}}$ intensities of primary (solid) and secondary (dashed) scattering at ${\boldsymbol \rho}$ along the $y$ axis for spherical scatterers of different sizes. For each simulation, we considered two scatterers of the same size spaced (${d_{\textit{qr}}}$) 4 µm apart. $d = 4\,\,\unicode{x00B5}{\rm m}$ and $\beta {= 45^ \circ}$. Four scatterer sizes considered are 0.5 µm, 1 µm, 1.5 µm, and 2 µm.
To obtain an estimate of the magnitude of the secondary scattering effect, we considered the scenario shown in Fig. 9(a) and computed the primary and secondary scattering at ${\boldsymbol \rho}$ along the $y$ axis for spherical scatterers. Figure 9(b) shows the primary and secondary scattering intensities for the arrangement in Fig. 9(a). In this scenario, we see that the secondary scattering intensities are, at most, two to three orders of magnitude smaller than the primary scattering intensities. Of course, these results do depend on the size and location of the secondary scatterer, but the scenario shown here is representative of the typical relative differences between primary and secondary scattering intensities.
In practice, the influence of tertiary or higher-order scattering can be safely ignored. We accordingly restrict our description to primary scattering and secondary scattering. With reference to Fig. 9(a), consider a field radiated from a dipole $j$ that is incident on a scattering object $Q$. In the case of a nearby object $R$, following Eq. (68) derived for a single scatterer, we can write the scattered field observed at point ${\boldsymbol \rho}$ for multiple scatterers as
Here, ${n_{{\rm scat}}}$ is the number of scatterers in the media, and ${\mathbb{S}_q}$ and ${\mathbb{S}_r}$ are the scattering amplitude matrices for the $Q$th and $R$th scatterers, respectively. ${\mathbb{R}_{2q\rho}}$ is a $2 \times 2$ transformation matrix at the $Q$th scatterer with an observation location at $\rho$. ${\mathbb{R}_{2r\rho}}$ is a $2 \times 2$ transformation matrix at the $R$th scatterer with an observation location at $\rho$. ${\mathbb{R}_{2qr}}$ is a $2 \times 2$ transformation matrix at the $Q$th scatterer with an observation location at the $R$th scatterer [Fig. 9(a)]. $[{\hat{\textbf{m}}},{\hat{\textbf{n}}}]_q^T$ and $[{\hat{\textbf{m}}},{\hat{\textbf{n}}}]_r^T$ transform parallel and perpendicular components to Cartesian components, and subscripts $q$ and $r$ represent $Q$th and $R$th scatterers, respectively.
When we consider all locations of the spherical reference surface, the total scattered field at ${\boldsymbol \rho}$ is given by
When the near-field Lorenz–Mie theory is considered, the $\exp (ikd)\mathbb{S}$ matrix in Eqs. (74) and (75) and the other primary and secondary scattering expressions in Part II [51] should be replaced with $\mathbb{E}(kd)$ as shown in Eq. (71).
Figure 10(a) shows the computed scattered electric field based on Eq. (76). Two scatterers are placed prior to the focal plane, and three scatterer diameters are considered: 1 µm, 2.5 µm, and 5 µm. The results are in excellent agreement with the FDTD results given in Fig. 10(b). We note that the results obtained from Eq. (76) may lack accuracy in the immediate near field of the scatterers because of the missing radial component. On the other hand, the FDTD results may be affected by the focused beam implementation technique, discretization errors, and unavoidable reflections from the perfectly matched layer (PML) [27,28]. The results obtained in Eq. (76) and FDTD simulations near the focal plane are compared numerically in [47].
Fig. 10. $x$ component of the scattered electric field amplitude on the $x - z$ plane ($y = 0$) computed using (a) our method discussed in Section 5 and (b) FDTD for 1.0 µm (left), 2.5 µm (middle), and 5.0 µm (right) spherical scatterers. The scatterers are located at (${-}{2.4}\;{\unicode{x00B5}{\rm m}}$, 0 µm, ${-}{6}\;{\unicode{x00B5}{\rm m}}$) and (2.4 µm, 0 µm, ${-}{15}\;{\unicode{x00B5}{\rm m}}$). The horizontal dashed line shows the focal plane. Amplitude results are shown in ${{\rm log}_{10}}$ scale. NA of the lens is 0.667.
When scatterers are involved, a detector at ${\boldsymbol \rho}$ captures the total electric field ${\textbf{ E}^{{\rm tot}}}({\boldsymbol \rho})$ representing the superposition of both incident $\textbf{ E}({\boldsymbol \rho})$ and scattered ${\textbf{ E}^s}({\boldsymbol \rho})$ fields [14,15,17]:
$\textbf{ E}({\boldsymbol \rho})$ in Eq. (77) is provided by Eq. (49) or Eq. (51). The corresponding intensity of the field at the location of a nonlinearly polarizable object is proportional to $|{E_x}{|^2} + |{E_y}{|^2} + |{E_z}{|^2}$, where ${E_x}$, ${E_y}$, and ${E_z}$ are Cartesian components of ${\textbf{ E}^{{\rm tot}}}({\boldsymbol \rho})$. It is the computation of the total electric field and corresponding intensities within the focal volume that provides the critical input to determine the generation of optical signals within the focal volume considered in Part II of this tutorial.
6. CONCLUSION AND FURTHER READING
We have provided a framework that utilizes classical analytical electromagnetic field propagation methods to comprehensively model optical microscopy in scattering samples with fixed scatterer configurations. The equations used in the tutorial for propagating light in non-scattering media originate from the vectorial and scalar theory of diffraction and focusing. Their foundation and application can be found in [15,39,40,53,54,57,64].
We have limited our scope to spherical scatterers, which are generally considered to provide a good approximation to the scattering objects relevant for optical microscopy of biological samples. Rigorous modeling and computation of the scattered field for non-spherical scatterers is more complicated and is treated in [48,84].
We have treated both primary and secondary scattering. The inclusion of secondary scattering increases the computational time exponentially. Secondary scattering becomes important once scatterers are larger than three to four times the wavelength and when they are located very close to each other. In other cases, the effect of secondary scattering on the final result is rather low and may be ignored. In cases where the scatterers are smaller than the wavelength with the neighboring scatterers located in the far zone, secondary scattering can be safely neglected.
In Part II of this tutorial [51], we consider focal field distributions in the presence of scatterers, signal generation, radiation, and far-field detection. We also provide a case study that includes two-photon and SHG microscopy.
Funding
National Institute of General Medical Sciences (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. B. A. Wilt, L. D. Burns, E. T. W. Ho, K. K. Ghosh, E. A. Mukamel, and M. J. Schnitzer, “Advances in light microscopy for neuroscience,” Annu. Rev. Neurosci. 32, 435–506 (2009). [CrossRef]
2. P. Kanchanawong and C. M. Waterman, “Advances in light-based imaging of three-dimensional cellular ultrastructure,” Curr. Opin. Cell Biol. 24, 125–133 (2012). [CrossRef]
3. C. Zhang, D. Zhang, and J.-X. Cheng, “Coherent Raman scattering microscopy in biology and medicine,” Annu. Rev. Biomed. Eng. 17, 415–445 (2015). [CrossRef]
4. V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods 7, 603–614 (2010). [CrossRef]
5. K. Wang, N. G. Horton, and C. Xu, “Going deep: brain imaging with multi-photon microscopy,” Opt. Photon. News 24(11), 32–39 (2013). [CrossRef]
6. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. 32, 2309–2311 (2007). [CrossRef]
7. D. Sinefeld, H. P. Paudel, D. G. Ouzounov, T. G. Bifano, and C. Xu, “Adaptive optics in multiphoton microscopy: comparison of two, three and four photon fluorescence,” Opt. Express 23, 31472–31483 (2015). [CrossRef]
8. K. M. Hampson, R. Turcotte, D. T. Miller, K. Kurokawa, J. R. Males, N. Ji, and M. J. Booth, “Adaptive optics for high-resolution imaging,” Nat Rev Methods Primers. 1, 68 (2021). [CrossRef]
9. O. Katz, E. Small, Y. Guan, and Y. Silberberg, “Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers,” Optica 1, 170–174 (2014). [CrossRef]
10. X. Tao, T. Lam, B. Zhu, Q. Li, M. R. Reinig, and J. Kubby, “Three-dimensional focusing through scattering media using conjugate adaptive optics with remote focusing (CAORF),” Opt. Express 25, 10368–10383 (2017). [CrossRef]
11. M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light Sci. Appl. 3, e165 (2014). [CrossRef]
12. J. Kubby, S. Gigan, and M. Cui, eds., Wavefront Shaping for Biomedical Imaging (Cambridge, 2019).
13. D. Débarre, M. J. Booth, and T. Wilson, “Image based adaptive optics through optimisation of low spatial frequencies,” Opt. Express 15, 8176–8190 (2007). [CrossRef]
14. H. C. van de Hulst, Light Scattering by Small Particles (Wiley, 1957).
15. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983).
16. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 1999).
17. G. Kristensson, Scattering of Electromagnetic Waves by Obstacles (SciTech, 2016).
18. 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]
19. 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, 043903 (2009). [CrossRef]
20. 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]
21. 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]
22. A. K. Glaser, Y. Chen, and J. T. C. Liu, “Fractal propagation method enables realistic optical microscopy simulations in biological tissues,” Optica 3, 861–869 (2016). [CrossRef]
23. 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]
24. 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]
25. 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]
26. 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]
27. 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]
28. I. R. Çapoglu, 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]
29. S. H. Tseng, “2-D PSTD Simulation of focusing monochromatic light through a macroscopic scattering medium via optical phase conjugation,” Biomed. Opt. Express 6, 815–826 (2015). [CrossRef]
30. D. Zhang, I. Capoglu, Y. Li, L. Cherkezyan, J. Chandler, G. Spicer, H. Subramanian, A. Taflove, and V. Backman, “Finite-difference time-domain-based optical microscopy simulation of dispersive media facilitates the development of optical imaging techniques,” J. Biomed. Opt. 21, 065004 (2016). [CrossRef]
31. I. R. Çapoglu, 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]
32. P. R. Munro, D. Engelke, and D. D. Sampson, “A compact source condition for modelling focused fields using the pseudospectral time-domain method,” Opt. Express 22, 5599–5613 (2014). [CrossRef]
33. Z. Wu, Y. Han, J. Wang, and Z. Cui, “Generation of Bessel beam sources in FDTD,” Opt. Express 26, 28727–28737 (2018). [CrossRef]
34. A. Mazzolani, C. M. Macdonald, and P. R. T. Munro, “Application of a Taylor series approximation to the Debye-Wolf integral in time-domain numerical electromagnetic simulations,” J. Opt. Soc. Am. A 39, 927–935 (2022). [CrossRef]
35. 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]
36. 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]
37. 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) (2017), pp. 1106–1108.
38. A. S. van de Nes, L. Billy, S. F. Pereira, and J. J. M. Braat, “Calculation of the vectorial field distribution in a stratified focal region of a high numerical aperture imaging system,” Opt. Express 12, 1281–1293 (2004). [CrossRef]
39. P. R. Munro, “Tool for simulating the focusing of arbitrary vector beams in free-space and stratified media,” J. Biomed. Opt. 23, 090801 (2018). [CrossRef]
40. M. R. Foreman and P. Török, “Computational methods in vectorial imaging,” J. Mod. Opt. 58, 339–364 (2011). [CrossRef]
41. J. P. Barton, D. R. Alexander, and S. A. Schaub, “Internal and near-surface electromagnetic fields for a spherical particle irradiated by a focused laser beam,” J. Appl. Phys. 64, 1632–1639 (1988). [CrossRef]
42. A. S. van de Nes and P. Torok, “Rigorous analysis of spheres in Gauss-Laguerre beams,” Opt. Express 15, 13360–13374 (2007). [CrossRef]
43. T. X. Hoang, X. Chen, and C. J. R. Sheppard, “Interpretation of the scattering mechanism for particles in a focused beam,” Phys. Rev. A 86, 033817 (2012). [CrossRef]
44. P. B. Bareil and Y. Sheng, “Modeling highly focused laser beam in optical tweezers with the vector Gaussian beam in the T-matrix method,” J. Opt. Soc. Am. A 30, 1–6 (2013). [CrossRef]
45. A. D. Kiselev and D. O. Plutenko, “Mie scattering of Laguerre-Gaussian beams: photonic nanojets and near-field optical vortices,” Phys. Rev. A 89, 043803 (2014). [CrossRef]
46. 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]
47. 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]
48. M. I. Mishchenko, L. D. Davis, and A. A. Lacis, Multiple Scattering of Light by Particles (Cambridge University, 2006).
49. D. W. Mackowski and M. I. Mishchenko, “A multiple sphere T-matrix Fortran code for use on parallel computer clusters,” J. Quant. Spectrosc. Radiat. Transfer 112, 2182–2192 (2011). [CrossRef]
50. 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]
51. J. C. Ranasinghesagara, E. O. Potma, and V. Venugopalan, “Modeling nonlinear optical microscopy in scattering media, part II. Radiation from focal volume to far-field: tutorial,” J. Opt. Soc. Am. A 40, 883–897 (2023). [CrossRef]
52. E. Hecht, Optics, 4th ed. (Addison-Wesley, 2001).
53. J. Braat and P. Török, Imaging Optics (Cambridge University, 2019).
54. P. Török, P. D. Higdon, and T. Wilson, “On the general properties of polarised light conventional and confocal microscopes,” Opt. Commun. 148, 300–315 (1998). [CrossRef]
55. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2006).
56. R. Azzam and N. Bashra, Ellipsometry and Polarized Light (Elsevier, 1987).
57. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems. II. Structure of the image field in an aplanatic system,” Proc. Royal Soc. London A 253, 358–379 (1959). [CrossRef]
58. J. Kim, Y. Wang, and X. Zhang, “Calculation of vectorial diffraction in optical systems,” J. Opt. Soc. Am. A 35, 526–535 (2018). [CrossRef]
59. J. J. Stamnes, Waves in Focal Regions: Propagation, Diffraction and Focusing of Light, Sound and Water Waves (Hilger, 1986).
60. H. H. Hopkins, “Herschel’s condition,” Proc. Phys. Soc. 58, 100–105 (1946). [CrossRef]
61. E. Abbe, “On the estimation of aperture in the microscope,” J. R. Microsc. Soc. 1, 388–423 (1881). [CrossRef]
62. V. S. Ignatowsky, “Diffraction by a lens of arbitrary aperture,” Trans. Opt. Inst. Petrograd 1, 1–36 (1919).
63. R. K. Luneburg, Mathematical Theory of Optics (University of California, 1966).
64. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).
65. J. A. Stratton and L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. 56, 99–107 (1939). [CrossRef]
66. P. Török, “Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number,” J. Opt. Soc. Am. A 15, 3009–3015 (2006). [CrossRef]
67. H. Piaggio, “The mathematical theory of Huygens’ principle,” Nature 145, 531–532 (1940). [CrossRef]
68. D. A. B. Miller, “Huygen’s wave propagation principle corrected,” Opt. Lett. 16, 1370–1372 (1991). [CrossRef]
69. F. Depasse, M. A. Paesler, D. Courjon, and J. M. Vigoureux, “Huygens-Fresnel principle in the near field,” Opt. Lett. 20, 234–236 (1995). [CrossRef]
70. E. Wolf and Y. Li, “Conditions for the validity of the Debye integral representation of focused fields,” Opt. Commun. 39, 205–210 (1981). [CrossRef]
71. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds., NIST Handbook of Mathematical Functions (Cambridge University, 2010).
72. C. G. Koay, “A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere,” J. Comput. Sci. 2, 377–381 (2011). [CrossRef]
73. M. O’Neil, “A generalized Debye source approach to electromagnetic scattering in layered media,” J. Math. Phys. 55, 012901 (2014). [CrossRef]
74. M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus field calculations,” Opt. Express 14, 11277–11291 (2006). [CrossRef]
75. M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express 12, 6530–6539 (2004). [CrossRef]
76. M. I. Mishchenko, “Far-field approximation in electromagnetic scattering,” J. Quant. Spectrosc. Radiat. Transfer 100, 268–276 (2006). [CrossRef]
77. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, 1941).
78. G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists, 7th ed. (Academic, 2012).
79. S. Boyd and L. Vandenberghe, Introduction to Linear Algebra: Vectors, Matrices and Least Squares (Cambridge University, 2018).
80. Z. J. Smith and A. J. Berger, “Validation of an integrated Raman- and angular-scattering microscopy system on heterogeneous bead mixtures and single human immune cells,” Appl. Opt. 48, D109–D120 (2009). [CrossRef]
81. M. Xu, T. T. Wu, and J. Y. Qu, “Unified Mie and fractal scattering by cells and experimental study on application in optical characterization of cellular and subcellular structures,” J. Biomed. Opt. 13, 024015 (2014). [CrossRef]
82. A. C. Lesina, J. van der Kolk, P. Berini, and L. Ramunno, “Computational electrodynamics - a powerful tool for nanophotonics and microscopy,” MRS Adv. 3, 753–760 (2018). [CrossRef]
83. F. Slimani, G. Grehan, G. Gouesbet, and D. Allano, “Near-field Lorenz-Mie theory and its application to microholography,” Appl. Opt. 23, 4140–4148 (1984). [CrossRef]
84. M. I. Mishchenko, “Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation,” Appl. Opt. 39, 1026–1031 (2000). [CrossRef]
