PyMieSim: an open-source library for fast and flexible far-field Mie scattering simulations

Martin Poinsinet de Sivry-Houle; Nicolas Godbout; Nicolas Godbout; Caroline Boudoux; Caroline Boudoux; Caroline Boudoux

doi:10.1364/OPTCON.473102

1. Introduction

Light-particle interactions are at the heart of a wide range of technologies, from biomedical imaging to solid-state physics [1,2]. In particular, numerous biomedical optics applications rely upon light scattering interactions for early diagnosis of cancer [3] or a variety of diseases, such as Alzheimer’s [4]. Understanding such complex interactions between light sources, scatterers, and detectors is of prime importance in validating experimental results and theoretical predictions. Optical configurations composed of a source, scatterer, and detector (SSD) are defined by several attributes resulting in systems of high dimensionality. Although the source and scatterer interaction can be challenging to solve, under certain conditions [5], the scatterer can be approximated to a geometry sporting symmetries such as a sphere, simplifying the mathematical expressions and, thus, enabling the use of fast simulation tools, such as the Lorenz-Mie scattering theory (LMT). The LMT framework allows finding exact solutions to Maxwell’s equations for some scatterer geometries such as homogeneous spheres, stratified spheres, or circular cylinders, albeit using planewave illumination. Bohren and Huffman’s Fortran codes [1] for such a framework have been widely used and adapted in several programming languages, providing an essential tool for experimental validation and theoretical predictions. Python is currently a prevalent, interpreted programming language [6] that has become pivotal for scientific applications, notably due to its wide range of available libraries [7,8], exportability, and flexibility. Presently, a few Python/C++ projects exist for LMT such as PyMieScatt (Python) [9], MieLab (Python) [10] or Bim-Sim (C++) [11]. These do not emulate interactions with detectors (coherent or incoherent), and nor are they suitable for studying optical properties of large sets of configurations. Moreover, as these projects are written exclusively in Python or C++, they inherit either the slowness or complexity of the associated programming language. In this paper, we present PyMieSim, a fast and flexible Python library wrapping a C++/Fortran core for simulating a wide range of SSD configurations with a much-reduced run time. The library allows defining a complete SSD system to compute the full scattered field (over 4 $\pi$ steradian), properties such as efficiencies (scattering, extinction, absorption, and radiation pressure), asymmetry, and albedo parameters, as well the coupling between the far-field scattered wave and a detection scheme (incoherent detector or coherent optical fibers). PyMieSim comprises three different LMT solvers, each one computing solutions of Maxwell’s equations for different scatterer form factors: homogeneous spheres, 2-layer spheres, and perpendicularly illuminated infinite circular cylinders. Whilst there exists a few libraries allowing for the study of source-scatterer interactions [9,12,13], none allow evaluating interaction with a detector. More importantly, other packages are not suitable for studying the optical properties of a large set of optical configurations with the required short compute time. By contrast, PyMieSim allows the simulation of coherent detection schemes like fiber optic modes or incoherent schemes like photodiode detectors. PyMieSim offers more than two orders of magnitude in reduced compute-time over all comparable publicly-available packages. In the following sections, we discuss the most important aspects of the scattering algorithm implemented in PyMieSim, validate its results against state-of-the-art algorithms, and demonstrate typical simulations, including various detection geometries.

1.1 Lorenz-Mie theory applications

Lorenz-Mie scattering theory (LMT) is a central framework for describing various processes that are at the core of several technologies and studies. LMT makes use of electric and magnetic multipole expansions to compute the scattered near- and far fields. Originally, LMT packages only allowed computing solutions for homogeneous spheres [14]. Since then, various geometries were adapted, including: infinite cylinders [15], stratified spheres [10], and others [16]. The analysis of such solutions as a function of the scattering event parameters may serve many applications in biomedical optics, such as light scattering spectroscopy (LSS) [17] or surface plasmon resonance enhancement, respectively. Solutions to Mie scattering may also be used for particle sizing [18], absorption quantification through the atmosphere [19], sensing of cells in biological tissues [20], flow cytometry [21], and also in metamaterial design [22]. While the LMT framework is bound to assume planewave illumination and is limited to certain scatterer geometries, it is widely used as a first-order approximation for more complex configurations. The generalized Lorenz-Mie scattering theory (GLMT) framework is an extension of LMT that enables defining an arbitrary illumination source to compute the scattering attribute. However, GLMT is more computationally expensive—especially when rotational symmetry is broken—as it usually involves computing higher-order terms of beam decomposition [23].

2. Theory

PyMieSim is a Python package that was developed to provide a novel tool for the extensive study of light scattering interactions. In order to simulate a scattering event as shown in Fig. 1, PyMieSim exploits the LMT framework. From this, the library can output various scattering properties, compute the scattered far fields and compute the coupling between the latter and different detection schemes.

Fig. 1. Representation of a scattering event configuration. An incident planewave (wavelength $\lambda _0$, polarization and electric amplitude $\vec E_0$, and wavevector $\hat k$) interacts with a homogeneous spherical scatterer (diameter $d$, and index, ${m}_{\textrm {sca}}$) immersed in a medium (index ${m}_{\textrm {med}}$). The resulting scattered wave (not shown) is projected onto an incoherent detector (numerical aperture NA, and angular location $\gamma$ and $\phi$) through a polarizing filter (not shown). Alternatively, PyMieSim allows replacing the incoherent detector with a coherent detection scheme (not shown) such as a single- or few-mode optical fiber (defined as a superposition of the scattered field onto several (LP) propagation modes).

Download Full Size | PDF

2.1 Optical setup

Figure 1 shows a typical geometry for a scattering event defined using variables found in the package. Illumination is described using the wavelength in vacuum, $\lambda _0$, the polarization, $\hat E_0$ and amplitude $\left | E_0 \right |$. A scatterer (either a homogeneous sphere, a stratified sphere, or an infinite cylinder) is described through its diameter, $d$, its refractive index, ${m}_{\textrm {sca}}$, and that of the surrounding medium, ${m}_{\textrm {med}}$. A detector defined with a (NA) is placed at an inclination and azimuthal angles $\phi$ and $\gamma$, respectively, accepts two orthogonally polarized scattered fields that can be specifically selected with a polarization filter. The detector coupling mechanism can also be set to coherent or incoherent, such as a photodiode or linearly polarized (LP) propagation modes in an optical fiber, respectively.

2.2 Lorenz-Mie theory

The standard LMT provides a mathematical framework that enables numerical methods to infer scattering properties of a scatterer planewave interaction. While approximate algorithms such as the anomalous diffraction [24] or Wiscombe approximation [25] were developed to accelerate the computation, the need for such speed up has been decreasing with the evolution of central process unit (CPU) and graphical process unit (GPU) architectures and clock speed. As of today, a well-written exact solution algorithm can usually achieve the time-restrictive requirements of the end-user. For this reason, we will restrain ourselves to discussing exact solution algorithms. In particular, LMT returns the scattering matrix elements, denoted as $S_{1-4}$, and defined as

(1)$$\begin{pmatrix}E^s_{{\parallel}} \\ E^s_{{\perp}} \end{pmatrix} = \frac{e^{ikr}}{ikr} . \begin{pmatrix}S_1 & S_4 \\ S_3 & S_2 \end{pmatrix} \begin{pmatrix}E^i_{{\parallel}} \\ E^i_{{\perp}} \end{pmatrix} \quad,$$

where $E^i$ and $E^s$ are the incident and scattered electric fields, respectively, $k$ is the wavenumber, and $r$ is the distance from the scatterer. For a spherical scatterer, the $S_3$ and $S_4$ are null. Other useful coefficients are defined in

(2)$$\vec{E}^s = E_0 \sum^{\infty}_{n=1} \frac{i^n (2n+1)}{n (n+1)} \left( i a_n \vec{N}^{(3)}_{e1n} - b_n \vec{M}^{(3)}_{o1n} \right) \quad,$$

where $\vec {N}^{(3)}_{e1n}$, and $\vec {M}^{(3)}_{o1n}$ are the vector spherical harmonics (VSH) as defined in [1]. LMT’s exact algorithms usually represent the scattered far-field amplitudes as an infinite summation of VSH associated with electric and magnetic multipole coefficients, $a_n$ and $b_n$, respectively. The coefficients $a_n$ and $b_n$ are used to compute scattering and absorption efficiencies, for example. The VSH, also denoted as electric and magnetic multipoles, form a complete orthogonal basis for the decomposition of the scattered fields. The scattering near-fields are also defined within the LMT and are associated to the coefficients $c_n$ and $d_n$, which are not of interest in this paper. Most of the computational time for those algorithms is used to evaluate the scattering coefficient up to a certain order $n_{max}$. The scattering coefficients ($a_n$ and $b_n$) are derived using the Legendre and Bessel functions of order $n$, which are computed recursively. For calculations that do not involve near-field computing, Wiscombe’s criterion is usually adopted [26]. This criterion defines the maximum order as

(3)$$n_{max} = x + 4.05 x^{1/3} + 2 \quad ,$$

where $x = \pi d / \lambda$ is the size parameter, $d$ is the diameter of the scatterer, and $\lambda$ is the incident wavelength in the surrounding medium ($\lambda =\lambda _0/{m}_{\textrm {med}}$). The homogeneous sphere, multilayer sphere and infinite cylinder are specific cases for which the scattering coefficients are well-defined as they lead to an analytical expression. In this paper, we follow the same mathematical notation as developed in [27]:

(4)$$\begin{aligned} \begin{cases} \Psi_n(x) & = x \sqrt{\frac{\pi}{2 x}} J_{n+1/2} (x) \\ \xi_n(x) & = x\sqrt{\frac{\pi}{2 x}} H^{(1)}_{n+1/2} (x) \end{cases}, \end{aligned}$$

with $J_n$, $H_n^{(1)}$ being the cylindrical Bessel and Hankel functions of the first kind, respectively. When a scatterer is also absorbing, its refractive index (RI) become complex, and so is the argument $x$ in Eq. (4). However, to the best of our knowledge, no C++ library implements complex argument Bessel functions. We therefore opted for a Fortran wrapper which exploits the Amos library [28] capable of such computations. Adopting this convention, and following Bohren and Huffman [1], we write the scattering coefficients $a_n$ and $b_n$ for the homogeneous sphere, infinite cylinder and 2-layer sphere as

(5)$$\textrm{homogeneous sphere} : \begin{cases} a_n = \frac{ \tilde{m} \Psi_n(\tilde{m} x) \dot{\Psi}_n (x) \,-\, \Psi_n (x) \dot{\Psi}_n(\tilde{m} x)} { \tilde{m} \Psi_n(\tilde{m} x) \dot{\xi}_n (x) \,-\, \xi_n (x) \dot{\Psi}_n(\tilde{m} x) } \\ b_n = \frac{ \Psi_n(\tilde{m} x) \dot{\Psi}_n (x) \,-\, \tilde{m} \Psi_n (x) \dot{\Psi}_n(\tilde{m} x)} { \Psi_n(\tilde{m} x) \dot{\xi}_n (x) \,-\, \tilde{m} \xi_n (x) \dot{\Psi}_n(\tilde{m} x) } \\ \end{cases},$$

(6)$$\textrm{infinite cylinder} : \begin{cases} a_n = \frac{ \tilde{m} \dot{J}_n(x) J_n (\tilde{m} x) \,-\, J_n(x) \dot{J}_n (\tilde{m} x) } { \tilde{m} J_n(\tilde{m} x) \dot{H}^{(1)}_n (x) \,-\, \dot{J}_n(\tilde{m} x) H^{(1)}_n (x) }\\ b_n = \frac{ \tilde{m} J_n(\tilde{m} x) \dot{J}_n (x) \,-\, \tilde{m} \dot{J}_n(\tilde{m} x) J_n (x) } { \tilde{m} J_n(\tilde{m} x) \dot{H}^{(1)}_n (x) \,-\, \tilde{m} \dot{J}_n(\tilde{m} x) H^{(1)}_n (x) } \end{cases},$$

(7)$$\textrm{2-layer sphere} : \begin{cases} a_n = \frac{ \Psi_n (y) \left[ \dot{\Psi}_n (\tilde{m}_{2} y) - A_n (x) \dot{\xi}_n (\tilde{m}_{2} y) \right] - \tilde{m}_{2} \dot{\Psi}_n (y) \left[ \Psi_n (\tilde{m}_{2} y) - A_n (x) \xi_n (\tilde{m}_{2} y) \right] } { \xi_n (y) \left[ \dot{\Psi}_n (\tilde{m}_{2} y) - A_n (x) \dot{\xi}_n (\tilde{m}_{2} y) \right] - \tilde{m}_{2} \dot{\xi}_n (y) \left[ \Psi_n (\tilde{m}_{2} y) - A_n (x) \xi_n (\tilde{m}_{2} y) \right] }\\ b_n = \frac{ \tilde{m}_{2} \Psi_n (y) \left[ \dot{\Psi}_n (\tilde{m}_{2} y) - B_n (x) \dot{\xi}_n (\tilde{m}_{2} y) \right] - \dot{\Psi}_n (y) \left[ \Psi_n (\tilde{m}_{2} y) - B_n (x) \xi_n (\tilde{m}_{2} y) \right] } { \tilde{m}_{2} \xi_n (y) \left[ \dot{\Psi}_n (\tilde{m}_{2} y) - A_n (x) \dot{\xi}_n (\tilde{m}_{2} y) \right] - \dot{\xi}_n (y) \left[ \Psi_n (\tilde{m}_{2} y) - B_n (x) \xi_n (\tilde{m}_{2} y) \right] } \end{cases},$$

\textrm{with} : \begin{cases} A_n = \frac{ \tilde{m}_{2} \Psi_n (\tilde{m}_{2} x) \dot{\Psi}_n (\tilde{m}_{1} x) - \tilde{m}_{1} \dot{\Psi}_n (\tilde{m}_{2} x) \Psi_n (\tilde{m}_{1} x)} {\tilde{m}_{2} \xi_n (\tilde{m}_{2} x) \dot{\Psi}_n (\tilde{m}_{1} x) - \tilde{m}_{1} \dot{\xi}_n (\tilde{m}_{2} x) \Psi_n (\tilde{m}_{1} x)}\\ B_n = \frac{\tilde{m}_{2} \psi_n (\tilde{m}_{1} x) \dot{\psi}_n (\tilde{m}_{2} x) - \tilde{m}_{1} \psi_n (\tilde{m}_{2} x) \dot{\psi}_n (\tilde{m}_{1} x)} {\tilde{m}_{2} \psi_n (\tilde{m}_{1} x) \dot{\xi}_n (\tilde{m}_{2} x) - \tilde{m}_{1} \xi_n (\tilde{m}_{2} x) \dot{\psi}_n (\tilde{m}_{1} x)} \end{cases} ,

where $x = k d_{\textrm {1}}/2$ and $y = k d_{\textrm {2}}/2$ are the size parameters for core and shell, respectively, with $d_{1,2}$ being the diameters of the associated scatterer layers, and $k$ the wavenumber. The dot ( $\dot {}$ ) represents a total derivative. The variable $\tilde {m}$ represents the relative RI of the scatterer, meaning that its value is divided by the RI of the medium. In the case of the 2-layer sphere, $\tilde {m}_1$ and $\tilde {m}_2$ represent the relative RI of the core and shell. Equations (5–7) assume a non-magnetic scatterer immersed in a non-absorbing medium.

2.3 Scattering properties

While the scattering coefficients represent the amplitudes for the VSH, which do not describe measurable properties of the scattering event, they can be used to compute the scattering efficiencies and other parameters. Table 1 shows PyMieSim’s output for each scatterer type. For the homogeneous and 2-layer sphere, PyMieSim computes the efficiencies for scattering ($Q_{\textrm {sca}}$), extinction ($Q_{\textrm {ext}}$), absorption ($Q_{\textrm {abs}}$), radiation pressure ($Q_{\textrm {pr}}$) and backscatter ($Q_{\textrm {back}}$), as well as the asymmetry parameter ($g$) and albedo ($\Lambda$), using the relations shown in the Table 1. These relations are based on the coefficients $a_n$ and $b_n$ described in Eqs. (5–7). However, due to the symmetry breaking for the infinite cylinder, the asymmetry parameter $g$ cannot be readily written in terms of the scattering coefficient $a_n$ and $b_n$. As such, PyMieSim evaluates the asymmetry parameter $g$ by computing the expected value of the inclination angle: $\left\langle\cos (\phi )\right\rangle$ (Equations not shown [1]).

Table 1. Scattering parameters computed by PyMieSim for spherical particles. PyMieSim also computes scattering parameters for infinite cylinders using equations found in [1] and reproduced in Code 1 [29], but not included here for conciseness’ sake.

View Table | View all tables in this article

The cross-sections, $\sigma$, associated with the efficiencies can be retrieved, by multiplying $Q$ with the geometrical cross-section of the scatterer $A$ as

\sigma_{i} = Q_i * A \quad [\textrm{m}^2].

2.4 Coupling mechanism

In addition to computing the scattering properties and scattering far-field, the PyMieSim package can evaluate the coupling of the scattering far-field with a specific detector. To do so, the solid angle space, where the far-field is evaluated, is generated using the Fibonacci lattice algorithm [30], which exhibits useful integration properties such as avoiding the local anisotropy of a conventional latitude–longitude array. The simulated detectors can be divided into two categories, namely, incoherent and coherent detection schemes. The first category embodies a system for which the output signal is proportional to the intensity of the impinging light field, such as a PIN diode. Such a detection scheme does not depend on the impinging field phase. The second category represents a system sensitive to the amplitude of the incident light field, such as a single-mode fiber (SMF) or a few-mode fiber (FMF). Furthermore, the coherent coupling mechanism can be subdivided into two parts: centered coupling and mean coupling. The first case occurs when the scatterer is perfectly centered with the detector in comparison to the incident light field. In that case, following the same argument described in [31], we define the centered coupling as

(8)$$\eta_{\textrm{centered}} = \left|{ \displaystyle \iint_{\Omega} \widetilde{\psi}(\xi,\nu)^* \, \widetilde{\phi}_{\textrm{det.}}(\xi,\nu) \, d\xi \, d\nu}\right|^2,$$

where $\widetilde {\psi }$ and $\widetilde {\phi }$ are the far-field representations of the scattered light field and the detector, respectively, $\xi$ and $\nu$ are angles in spherical coordinates, and, finally, $\Omega$ is the solid angle subtended by the NA of the imaging system. The second case occurs when the scatterer is allowed not to be centered (said to be offsetted) with the detector, we thus define the mean coupling.

(9)$$\eta_{\textrm{mean}} = \Big\langle \,\, \Big|\, \widetilde{\psi}(\xi,\nu) \, . \, \widetilde{\phi}_{\textrm{det.}}(\xi,\nu) \Big|^2 \, \Big\rangle_\Omega ,$$

For incoherent coupling mechanisms, the centered and mean couplings lead to the same results, as it does not depend on the phase of the scattering far field. Additionally, the PyMieSim software can also compute the scatterer footprint, which represents the coupling as a function of the scatterer spatial offset in two directions. Such a footprint is what we can expect to see using a scanning system, as it was demonstrated in [31].

3. Main features of PyMieSim

PyMieSim was developed using an architecture promoting fast computing time and ease of operation under several platforms.

3.1 Architecture

Figure 2 represents the structure of PyMieSim. The package takes advantage of the object-oriented aspect of Python to define three C++ classes: Source describes the incoming wave, Scatterer defines the scatterer geometry and attributes, and Detector defines the detector type, position, and attributes. Within the Python wrapper Experiment, Source and Scatterer interact together to return scattering efficiencies, while Detector returns the coupling ratio and scatterer footprint. The wrapping of these classes into a Python object is done using the PyBind11 library [32]. PyMieSim was developed in parallel with another Python package: PyOptik [33] to define the scatterer and immersion medium as specific materials with wavelength-dependent properties fetched from an online database (RefractiveIndex.info [34]). The code being composed of a Python wrapper for C++, it has to be compiled and, as such, is operating system (OS) dependent. For portability, PyMieSim exploits the Python package manager with uploaded pre-compiled code (wheel format). PyMieSim is currently available on the three main OS: Linux, macOS, Windows for different Python versions: 3.7, 3.8, 3.9. Furthermore, the project also allows end-users to manually install the package for any system configuration. PyMieSim is also available on the online Google Colab platform and, thus, may be used on any OS.

Fig. 2. Representation of the structure and computing capabilities of the PyMieSim library, highlighting the numerous attributes defining the optical configurations.

Download Full Size | PDF

3.2 Compute time

Wiscombe’s criterion, used to define the coefficient order, $n$, (see Eq. (3)) shows that the compute time rapidly increases for large size parameters, $x$. PyMieSim, written in C++ and Fortran, significantly reduces the runtime compared to existing Python packages, allowing for the study of a large set of SSD configurations up to large-size parameters. Figure 3(a) shows the runtime comparison for scattering coefficient calculations using different software packages. The calculations were performed on a Linux distribution using an AMD Ryzen 5 2600 (6 cores) CPU and 16 GB of ddr4 ram, and show a gain exceeding 130x for large sets of parameters (> 40 000). PyMieSim’s short computational time allows interfacing with other algorithms for fast post-processing, such as Mie inverse problems [35,36]. Such gain permits optimization tools to be used in conjunction with the library.

Fig. 3. Library bench-marking presenting the computing time for scattering coefficient $An$ for scatterer with RI of 4, size ranging from 400 to 1000 nm and a source with a wavelength of 800 nm.

Download Full Size | PDF

3.3 User interface

PyMieSim simulates SSD configurations, which are composed of up to fourteen attributes, as shown in Fig. 2. The package gives to the end-user full control over each of those attributes, allowing to iterate over ranges for each one. Moreover, PyMieSim comes with an additional tool (DataVisual [37]) for high-dimensionality data operation, analysis, and visualization. PyMieSim provides a handy application programming interface (API), accompanied by visual tools at each level to grasp the optical configurations. PyMieSim is also a set of independent tools that allow the user to embed within another project. The software implements a seamless integration of the C++ outputted data with the scientific NumPy library [38], allowing the package to be readily implemented in any other Python project. Figure 4 shows the simulation of the scattering event as depicted in Fig. 1, for ranging values of detector’s inclination angle $\phi$ and numerical aperture NA. The plots show that, as the NA of the detector increases, finer details in the scattering diagram vanish. Normalization is necessary to visualize all four curves, as increased NA results in increased coupling power, which the package allows calculating (not shown). Working examples of the code used to generate the plots shown in this paper can also be found at Code 1 [29] and https://pymiesim.readthedocs.io/en/latest/. As such, simulation and visualization are important aspects for an exhaustive understanding of the involved mechanisms. PyMieSim makes use of different graphical tools to give better insight to the end user. Fig. 5 shows examples of simulations performed using PyMieSim. Figure 5(a)–(c) show scattering experiments using a single scatterer ($d = {1000}\;\textrm{nm}$, $m=1.4$) and a linearly polarized planewave ($\lambda _0 = {500}\;\textrm{nm}$). Figure 5(a) shows amplitude scattering matrix elements $S_1(\phi )$ (left) and $S_2(\phi )$ (right). Figure 5(b) shows a 3D rendering of the real and imaginary parts of the scattering phase function, while Fig. 5(c) shows the far-field representation of the scattered intensity. Figure 5(d) shows different detector far-field representations. From left to right: an incoherent detector (NA=0.5), and three fiber modes: LP$_{01}$, LP$_{11}$, and LP$_{21}$. The LP far-field modes are computed using the Fraunhofer approximation.

Fig. 4. Results from the scattering experiment presented in Fig. 1 with $\lambda _0=1200$nm, $d=2000$nm, and a particle composed of BK7 glass. Coupling power is evaluated for four different detector NA’s displaced along the $\phi$ angle (and for $\gamma =0$).

Download Full Size | PDF

Fig. 5. Images obtained using the PyMieSim library for scattering simulations of various scatterer sizes illuminated with a 1000 nm planewave. (a) The non-null scattering matrix element S$_1$ and S$_2$; (b) Scattering phase function; (c & d) Scattering far-field in the $\hat {\phi }$ and $\hat {\gamma }$ direction, respectively. The real and imaginary parts are depicted in the right and left parts of the sub-figure; (e) Stokes parameters defining the polarization state using I (intensity), Q (linearly polarized, $0^\circ$), U (linearly polarized, $45^\circ$) and V (circularly polarized); (f) Detector far-field representation for a photodiode, LP$_{01}$, LP$_{11}$ and LP$_{21}$.

Download Full Size | PDF

4. Validation

Mie simulation results for the homogeneous sphere, the 2-layer sphere, and the infinite cylinder were extensively tested and compared to other libraries, namely: MiePlot [12], PyMieScatt [9], and MiePython [13]. However, to the best of our knowledge, no library implements the scattering calculation for infinite cylinders. As such, we compared our results to the literature. Figure 6 compares the scattering efficiencies for two orthogonal polarization, $Q_{scaI}, Q_{scaII}$, for a cylinder ($m=1.55$) in air (red and blue curves) and in water (orange and green curves) with results published in [1,39] (black dashed curves). The four curves overlap perfectly.

Fig. 6. Simulation of the scattering cross-sections for an infinite cylinder as a function of diameter. The simulation parameters are chosen to replicate the results presented in Fig.(8.7 & 8.8) of [1], which are here represented by the black dashed lines.

Download Full Size | PDF

5. Application

Here, we describe two use cases exploiting the unique features of PyMieSim, namely its fast compute core and the ability to simulate a coherent detection scheme.

5.1 Flow-cytometry optimization

PyMieSim’s fast computing core allows using the library for optimization problems. For example, we demonstrate how PyMieSim could be used to improve the accuracy of a flow cytometry system through the optimization of its detectors’ position. In its standard implementation, flow-cytometry uses two orthogonal photodetectors defined by the coordinates $\{ { \phi _0 = 0^\circ, \gamma _0=0^\circ }\}$ and $\{ { \phi _1 = 90^\circ, \gamma _1=0^\circ }\}$ in order to measure front and side scattering efficiencies of scatterers flowing one at a time through a monochromatic, polarized laser beam [40]. In this flow-cytometry optimization example, we define the optical system’s baseline parameters as in [21]. We thus modeled scattering from spherical scatterers within a phase space defined by a refractive index ranging from 1.33 to 1.6 and diameter ranging from 100 to 500 nm, illuminated with a polarized light source at 405 nm. We first modeled an unoptimized orthogonal detector pair with an NA of 0.15, then let PyMieSim optimize the detectors’ angular values ($\phi$ and $\gamma$) and NA (where NA is a free parameter ranging from [0.1 to 0.4]).

Figures 7(a-b) show, for each detector, the coupling power as a function of scatterer size and refractive index. For each detector, the measured coupling power defines an iso-contour, the intersection of which returns the scatterer’s size and refractive index. Uncertainty in the detector’s measurement broadens the iso-contour into areas bounded by dashed lines in Figs. 7(a)–(b) and shown as shaded areas in Figs. 7(c)–(d). Figures 7(c) maps the intersection of the areas in the unoptimized case. As iso-contours from each detector are not orthogonal, the uncertainty on size and refractive index remains large. PyMieSim was then used to compute a new set of detector’s angular locations and NA for which uncertainty would be minimized across all scatterer sizes and refractive indices. For a fixed measurement uncertainty, the best-case scenario occurs when the iso-contours for both detectors are orthogonal, this would mean that the projection of the gradient of the coupling field for each detector equals zero at every point. This was used as the optimization metric, normalized on a scale from zero to one, where zero means that all the iso-contour lines are orthogonal, and one means that they are co-linear. Mathematically, the optimization problem can be written in the following form:

(10)$$\min_{\forall s\, \in \, \boldsymbol{S}} \, \big(\vec{\nabla} F_0(s)\big) \cdot \big(\vec{\nabla} F_1(s)\big)\quad ,$$

where $\boldsymbol {S}$ is the feasible set $\{\phi \in \left [0, 360 \right ], \gamma \in \left [0, 360 \right ], \textrm {NA} \in \left [0.1, 0.4 \right ] \}$ as defined above, F$_0$ and F$_1$ are the coupling fields computed using PyMieSim for each detector. A quick assessment of the metric function establishes that the optimization problem is non-linear and non-convex, and thus we are not guaranteed to find a global minimum. Treating the question as a black-box optimization problem, we consequently choose the Nelder-Mead algorithm [41]. Figure 7(d) presents the iso-contours for the optimized detector’s angular values and NA, showing a reduction of uncertainty in size and refractive index. The optimization process gives us a local minimum of 0.27, while using the standard parameters, we obtain a metric value of 0.78. For such a demanding optimization process, PyMieSim’s short compute-time and easy-to-handle API were proven very useful for the reasonable computation time and expressiveness of a few lines of code.

Fig. 7. Parameter space optimization for a flow-cytometry experiment. (a-b) Coupling power of scattered light with a typical flow-cytometry system composed of two orthogonal detectors: (a) detector 0 $\{ { \phi _0 = 0^\circ, \, \gamma _0=0^\circ, \, \text {NA} = 0.15}\}$ (b) detector 1 $\{ {\phi _1 = 90^\circ, \, \gamma _1=0^\circ, \, \text {NA} = 0.15 }\}$. (c-d) Iso-contour projections of a specific scatterer measurement (RI: 1.48, diameter: 472 nm, dashed lines) with a measurement error of 5% (shaded areas) for the typical (c) and optimized (d) flow-cytometer, where optimized detectors are placed at: detector 0 $\{ {\phi _0 = 0^\circ, \, \gamma _0=42.3^\circ, \, \text {NA} = 0.38}\}$ and detector 1$\{{ \phi _1 = 133.2^\circ, \, \gamma _1=68.0^\circ, \, \text {NA} = 0.38 }\}$, respectively. $\Delta D$ and $\Delta n$ represent the diameter and refractive index uncertainty, respectively.

Download Full Size | PDF

5.2 Coherent imaging analysis

Another valuable application of the library concerns coherent or partially coherent imaging. PyMieSim allows simulating coherent detection schemes based on a fiber’s LP modes to improve endoscopy methods such as optical coherence tomography (OCT). Recently, the use of several LP modes has been under investigation as a way to acquire additional contrast [4,31,42]. PyMieSim allows modeling the source of contrast based on scatterer size and geometries for each LP mode [31]. Figure 8 shows an example of the coupling power of a coherent imaging system as a function of the scatterer size and different LP modes. In particular, Fig. 8 demonstrates that the coupling behavior does not change for small scatterers (i.e., approaching the Rayleigh regime), however, the coupling power in each mode diverges for larger scatterers, hinting at possible size-based contrast for FM-OCT.

Fig. 8. Simulating few-mode optical coherence tomography (FM-OCT) imaging at $\lambda _0=930nm$ and $\Delta \lambda =100nm$ with 3 LP modes. (a) Mean (solid lines) and standard deviation (shaded areas) of the power coupled back into each mode as a function of the scatterer diameter for LP01 (b), LP02 (c), and LP03 (d), represented in the far-field. Inset shows the overlap of the three coupling powers in the Rayleigh regime.

Download Full Size | PDF

6. Conclusion

We conclude by briefly comparing PyMieSim to other implementations of Mie theory. Table 2 summarizes the features of several Mie codes available online. The table highlights PyMieSim’s main aspects: enabling simulation for different types of scatterers, a fast computing core, an easy-to-install package, and an intuitive API. Moreover, PyMieSim addresses shortcomings associated with the other libraries, providing a continuous integration pipeline for an extensive testing process to ensure the library structure’s integrity and that the results produced match the literature. Additionally, the library provides complete online documentation and a platform with tens of examples featuring the reproduction of already published results.

Table 2. Features of open source LMT packages. Features were deduced from an in-depth analysis of each project.

View Table | View all tables in this article

The PyMieSim package is still under active development and has been designed to be readily upgraded. The library also makes use of external tools such as Travis [43] for continuous integration and Docker [44] for the compilation of the source code for each OS. There are various incremental additions and features that are being designed for near-future deployments, namely: GLMT framework, GPU computation, and graphical user interface (GUI). The GLMT framework is an extension of the LMT for structured illumination sources, thus allowing the user to define specific light sources such as a Gaussian beam [23], Bessel beam [45] or fiber optic LP modes. While many simulation projects exist for LMT simulations, only a few include GLMT. The framework can also be used to simulate the interaction with scatterers breaking the symmetries, as it was demonstrated for arbitrarily oriented spheroids [46]. While the necessary computation differs from the LMT framework, using the existing PyMieSim platform makes such a feature easier to implement. A similar approach could also be used to compute the near-field components of scattering events, allowing for multi-scatterer interaction and thus widening the range of applications of the library.

The ability to create experiment setups with ranging attributes for the scattering event is one of the main appeals of PyMieSim. Moreover, its implementation has been parallelized to make the most of the multiple cores that recent CPU architecture offers. Due to its highly independent implementation, the algorithm for such computation could be passed to GPU for an even more reduced runtime.

Finally, the aim of the package is to be easy to install and use. While the interpreted Python language allows using easy to handle platforms such as Jupyter notebook or Google Colab, the coding part might still be a barrier for some users. As such, developing a graphical interface can help overcome this difficulty. This is why we plan to develop a server-based application that will be accessible through a web browser.

Funding

Fonds de recherche du Québec – Nature et technologies (304141); Natural Sciences and Engineering Research Council of Canada (RGPIN-06404-2016, RGPIN-2018-06151).

Acknowledgements

We thank Pierre-Alexandre Aubé for his precious help with the writing of the PyMieSim package. Portions of this work were presented at SPIE Photonics West’s Symposium on Optical Interactions with Tissues and Cells XXXIII; and Advanced Photonics in Urology in 2022 (PC1195802) [47]. The presentation focused on the application programming interface and the usability of the package, while this manuscript delves into the equations behind the package, validates its output against other packages, demonstrates improvement in the computation time, and demonstrates the novel and unique features of PyMieSim through two use-cases of practical interest to the community.

Disclosures

The authors declare no conflicts of interest.

Data Availability

The source code for PyMieSim is available as Code 1 [29]. The raw data used in the article is available from the authors upon reasonable request.

References

1. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (John Wiley & Sons, 2008).

2. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, absorption, and emission of light by small particles (Cambridge university press, 2002).

3. S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” J. Biomed. Opt. 7(3), 329–340 (2002). [CrossRef]

4. P. Eugui, A. Lichtenegger, M. Augustin, D. J. Harper, M. Muck, T. Roetzer, A. Wartak, T. Konegger, G. Widhalm, C. K. Hitzenberger, A. Woehrer, and B. Baumann, “Beyond backscattering: optical neuroimaging by brad,” Biomed. Opt. Express 9(6), 2476–2494 (2018). [CrossRef]

5. M. Hammer, D. Schweitzer, B. Michel, E. Thamm, and A. Kolb, “Single scattering by red blood cells,” Appl. Opt. 37(31), 7410–7418 (1998). [CrossRef]

6. G. Van Rossum and F. L. Drake Jr, Python tutorial (Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995).

7. T. E. Oliphant, “Python for scientific computing,” Comput. Sci. Eng. 9(3), 10–20 (2007). [CrossRef]

8. K. J. Millman and M. Aivazis, “Python for scientists and engineers,” Comput. Sci. Eng. 13(2), 9–12 (2011). [CrossRef]

9. B. J. Sumlin, W. R. Heinson, and R. K. Chakrabarty, “Retrieving the aerosol complex refractive index using pymiescatt: A mie computational package with visualization capabilities,” J. Quant. Spectrosc. Radiat. Transf. 205, 127–134 (2018). [CrossRef]

10. O. Pe na-Rodríguez, P. P. G. Pérez, and U. Pal, “Mielab: a software tool to perform calculations on the scattering of electromagnetic waves by multilayered spheres,” Int. J. Spectrosc. 2011, 1–10 (2011). [CrossRef]

11. S. Berisha, T. van Dijk, R. Bhargava, P. S. Carney, and D. Mayerich, “Bim-sim: Interactive simulation of broadband imaging using mie theory,” Front. Phys. 5, 5 (2017). [CrossRef]

12. P. Laven, “MiePlot (a computer program for scattering of light from a sphere using Mie theory & the Debye series),” (2011), http://www.philiplaven.com/mieplot.htm.

13. S. Prahl, “MiePython,” (2021), https://miepython.readthedocs.io/en/latest/.

14. G. Mie, “Beiträge zur optik trüber medien, speziell kolloidaler metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]

15. W. Sun, N. G. Loeb, and B. Lin, “Light scattering by an infinite circular cylinder immersed in an absorbing medium,” Appl. Opt. 44(12), 2338–2342 (2005). [CrossRef]

16. F. Frezza, F. Mangini, and N. Tedeschi, “Introduction to electromagnetic scattering: tutorial,” J. Opt. Soc. Am. A 35(1), 163–173 (2018). [CrossRef]

17. V. Turzhitsky, L. Qiu, I. Itzkan, A. A. Novikov, M. S. Kotelev, M. Getmanskiy, V. A. Vinokurov, A. V. Muradov, and L. T. Perelman, “Spectroscopy of scattered light for the characterization of micro and nanoscale objects in biology and medicine,” Appl. Spectrosc. 68(2), 133–154 (2014). [CrossRef]

18. M. R. Vaziri, A. Omidvar, B. Jaleh, and N. P. Shabestari, “Investigating the extrinsic size effect of palladium and gold spherical nanoparticles,” Opt. Mater. 64, 413–420 (2017). [CrossRef]

19. T. C. Bond, G. Habib, and R. W. Bergstrom, “Limitations in the enhancement of visible light absorption due to mixing state,” J. Geophys. Res.: Atmos. 111(D20), D20211 (2006). [CrossRef]

20. P. Bassan, A. Kohler, H. Martens, J. Lee, E. Jackson, N. Lockyer, P. Dumas, M. Brown, N. Clarke, and P. Gardner, “Rmies-emsc correction for infrared spectra of biological cells: Extension using full mie theory and gpu computing,” J. Biophotonics 3(8-9), 609–620 (2010). [CrossRef]

21. E. van der Pol, L. de Rond, F. A. Coumans, E. L. Gool, A. N. Böing, A. Sturk, R. Nieuwland, and T. G. van Leeuwen, “Absolute sizing and label-free identification of extracellular vesicles by flow cytometry,” Nanomedicine (N. Y., NY, U. S.) 14(3), 801–810 (2018). [CrossRef]

22. Q. Zhao, J. Zhou, F. Zhang, and D. Lippens, “Mie resonance-based dielectric metamaterials,” Mater. Today 12(12), 60–69 (2009). [CrossRef]

23. J. Qiu and J. Shen, “Beam shape coefficient calculation for a gaussian beam: localized approximation, quadrature and angular spectrum decomposition methods,” Appl. Opt. 57(2), 302–313 (2018). [CrossRef]

24. H. C. Hulst and H. C. van de Hulst, Light scattering by small particles (Courier Corporation, 1981).

25. W. J. Wiscombe, “Improved mie scattering algorithms,” Appl. Opt. 19(9), 1505–1509 (1980). [CrossRef]

26. A. A. R. Neves and D. Pisignano, “Effect of finite terms on the truncation error of mie series,” Opt. Lett. 37(12), 2418–2420 (2012). [CrossRef]

27. G. Gouesbet and G. Gréhan, Generalized Lorenz-Mie theories, vol. 31 (Springer, 2011).

28. J. L. Arbuckle, Amos™ 7.0 user’s guide (Amos Development Corporation, 1995).

29. M. Poinsinet de Sivry-Houle and C. Boudoux, “MartinPdeS/PyMieSim: PyMieSim beta,” Zenodo (2022), https://doi.org/10.5281/zenodo.7407114.

30. R. Swinbank and R. James Purser, “Fibonacci grids: A novel approach to global modelling,” Q.J.R. Meteorol. Soc. 132(619), 1769–1793 (2006). [CrossRef]

31. M. Poinsinet de Sivry-Houle, S. B. Beaudoin, S. Brais-Brunet, M. Dehaes, N. Godbout, and C. Boudoux, “All-fiber few-mode optical coherence tomography using a modally-specific photonic lantern,” Biomed. Opt. Express 12(9), 5704–5719 (2021). [CrossRef]

32. W. Jakob, J. Rhinelander, and D. Moldovan, “pybind11–seamless operability between c++ 11 and python, GitHub (2017). https://github.com/pybind/pybind11

33. M. Poinsinet de Sivry-Houle and P.-A. Aubé, “PyOptik: A python package for material refractive index library,” Python Package Index (2021), https://pypi.org/project/PyOptik/.

34. M. N. Polyanskiy, “Refractive index database,”.

35. I. Ludlow and J. Everitt, “Inverse mie problem,” J. Opt. Soc. Am. A 17(12), 2229–2235 (2000). [CrossRef]

36. K. J. Chalut, M. G. Giacomelli, and A. Wax, “Application of mie theory to assess structure of spheroidal scattering in backscattering geometries,” J. Opt. Soc. Am. A 25(8), 1866–1874 (2008). [CrossRef]

37. M. Poinsinet de Sivry-Houle, “DataVisual: a python tool for high-dimensionality dataset visualization,” Python Package Index (2021), https://pypi.org/project/DataVisual/.

38. C. R. Harris, K. J. Millman, S. J. van der Walt, et al., “Array programming with NumPy,” Nature 585(7825), 357–362 (2020). [CrossRef]

39. T. A. Bashkatova, A. N. Bashkatov, V. I. Kochubey, and V. V. Tuchin, “Light-scattering properties for spherical and cylindrical particles: a simple approximation derived from mie calculations,” in Saratov Fall Meeting 2000: Optical Technologies in Biophysics and Medicine II, vol. 4241 (International Society for Optics and Photonics, 2001), pp. 247–259.

40. K. M. McKinnon, “Flow cytometry: an overview,” Curr. Protocols Immunology 120(1), 5 (2018). [CrossRef]

41. J. A. Nelder and R. Mead, “A simplex method for function minimization,” The Computer J. 7(4), 308–313 (1965). [CrossRef]

42. M. Kreysing, D. Ott, M. J. Schmidberger, O. Otto, M. Schürmann, E. Martín-Badosa, G. Whyte, and J. Guck, “Dynamic operation of optical fibres beyond the single-mode regime facilitates the orientation of biological cells,” Nat. Commun. 5(1), 5481 (2014). [CrossRef]

43. M. Beller, G. Gousios, and A. Zaidman, “TravisTorrent: synthesizing Travis CI and GitHub for full-stack research on continuous integration,” in 2017 IEEE/ACM 14th International Conference on Mining Software Repositories (MSR), (IEEE, 2017), pp. 447–450.

44. D. Merkel, “Docker: lightweight Linux containers for consistent development and deployment,” Linux J. 2014(239), 2 (2014). [CrossRef]

45. L. A. Ambrosio and H. E. Hernández-Figueroa, “Integral localized approximation description of ordinary Bessel beams and application to optical trapping forces,” Biomed. Opt. Express 2(7), 1893–1906 (2011). [CrossRef]

46. F. Xu, K. Ren, G. Gouesbet, G. Gréhan, and X. Cai, “Generalized Lorenz-Mie theory for an arbitrarily oriented, located, and shaped beam scattered by a homogeneous spheroid,” J. Opt. Soc. Am. A 24(1), 119–131 (2007). [CrossRef]

47. M. Poinsinet de Sivry-Houle, P.-A. Aube, N. Godbout, and C. Boudoux, “Pymiesim: an open source library for fast and flexible light scattering simulations,” in Optical Interactions with Tissue and Cells XXXIII; and Advanced Photonics in Urology, (SPIE, 2022), p. PC1195802.

	homogeneous sphere	2-layer sphere
$Q_{sca}$	$\frac{2}{x^{2}} \sum_{n} (2 n + 1) (\| a_{n} \|^{2} + \| b_{n} \|^{2})$
$Q_{ext}$	$\frac{2}{x^{2}} \sum_{n} (2 n + 1) R e {a_{n} + b_{n}}$
$Q_{abs}$	$Q_{e x t} - Q_{s c a}$
$Q_{back}$	$x^{- 2} \| \sum_{n} (2 n + 1) (- 1)^{n} (a_{n} + b_{n}) \|^{2}$
$Q_{ratio}$	$Q_{b a c k} / Q_{s c a}$
$Q_{pr}$	$Q_{e x t} - g Q_{s c a}$
$g = ⟨ \cos (ϕ) ⟩$	$\frac{4}{Q_{s c a} x^{2}} [\sum_{n} \frac{n (n + 2)}{(n + 1)} R e {a_{n} a_{n + 1}^{} + b_{n} b_{n + 1}^{}} + \sum_{n} \frac{(2 n + 1)}{n (n + 1)} R e {a_{n} b_{n}^{*}}]$
$Λ$	$Q_{e x t} / Q_{s c a}$

		PyMieSim	PyMieScatt	MieLab	MiePython	BimSim
	Homogeneous sphere	✓	✓	✓	✓	✓
	2-layer sphere	✓	✓	✓
	N-layer sphere			✓
Scatterer type	Infinite cylinder	✓
	Compiled core	✓			✓	✓
	CPU Parallel computing	✓
Performance	GPU computing					✓
	Python Interface	✓	✓	✓	✓
	Documentation	✓	✓	✓	✓
	Examples	✓	✓
	Package manager	✓	✓		✓
	Multi-platform	✓	✓		✓	✓
	Online platform	✓			✓
	Data visualization tools	✓
	Ranging parameters	✓
Ease of use	Graphic interface			✓
	Continuous Integration	✓
	Automatic testing	✓
Extra	Detector emulation	✓

	homogeneous sphere	2-layer sphere
$Q_{sca}$	$\frac{2}{x^{2}} \sum_{n} (2 n + 1) (\| a_{n} \|^{2} + \| b_{n} \|^{2})$
$Q_{ext}$	$\frac{2}{x^{2}} \sum_{n} (2 n + 1) R e {a_{n} + b_{n}}$
$Q_{abs}$	$Q_{e x t} - Q_{s c a}$
$Q_{back}$	$x^{- 2} \| \sum_{n} (2 n + 1) (- 1)^{n} (a_{n} + b_{n}) \|^{2}$
$Q_{ratio}$	$Q_{b a c k} / Q_{s c a}$
$Q_{pr}$	$Q_{e x t} - g Q_{s c a}$
$g = ⟨ \cos (ϕ) ⟩$	$\frac{4}{Q_{s c a} x^{2}} [\sum_{n} \frac{n (n + 2)}{(n + 1)} R e {a_{n} a_{n + 1}^{} + b_{n} b_{n + 1}^{}} + \sum_{n} \frac{(2 n + 1)}{n (n + 1)} R e {a_{n} b_{n}^{*}}]$
$Λ$	$Q_{e x t} / Q_{s c a}$

		PyMieSim	PyMieScatt	MieLab	MiePython	BimSim
	Homogeneous sphere	✓	✓	✓	✓	✓
	2-layer sphere	✓	✓	✓
	N-layer sphere			✓
Scatterer type	Infinite cylinder	✓
	Compiled core	✓			✓	✓
	CPU Parallel computing	✓
Performance	GPU computing					✓
	Python Interface	✓	✓	✓	✓
	Documentation	✓	✓	✓	✓
	Examples	✓	✓
	Package manager	✓	✓		✓
	Multi-platform	✓	✓		✓	✓
	Online platform	✓			✓
	Data visualization tools	✓
	Ranging parameters	✓
Ease of use	Graphic interface			✓
	Continuous Integration	✓
	Automatic testing	✓
Extra	Detector emulation	✓

PyMieSim: an open-source library for fast and flexible far-field Mie scattering simulations

Abstract

1. Introduction

1.1 Lorenz-Mie theory applications

2. Theory

2.1 Optical setup

2.2 Lorenz-Mie theory

2.3 Scattering properties

2.4 Coupling mechanism

3. Main features of PyMieSim

3.1 Architecture

3.2 Compute time

3.3 User interface

4. Validation

5. Application

5.1 Flow-cytometry optimization

5.2 Coherent imaging analysis

6. Conclusion

Funding

Acknowledgements

Disclosures

Data Availability

References

Supplementary Material (1)

Data Availability

Cited By

Figures (8)

Tables (2)

Equations (12)

Optics Continuum

Martin Poinsinet de Sivry-Houle	https://orcid.org/0000-0002-3637-7673
Nicolas Godbout	https://orcid.org/0000-0001-9494-1175
Caroline Boudoux	https://orcid.org/0000-0002-2145-9086