Toward Waveguide-Based Optical Chromatography

We report analytical expressions for optical forces acting on particles inside waveguides. The analysis builds on our previously reported Fourier Transform method to obtain Beam Shape Coefficients for any beam. Here we develop analytical expressions for the Beam Shape Coefficients in cylindrical and rectangular metallic waveguides. The theory is valid for particle radius a ranging from the Rayleigh regime to large microparticles, such as aerosols like virus loaded droplets. The theory is used to investigate how optical forces within hollow waveguides can be used to sort particles in “optical chromatography” experiments in which particles are optically propelled along a hollow-core waveguide. For Rayleigh particles, the axial force is found to scale with a⁶, while the radial force, which prevents particles from crashing into the waveguide walls, scales with a³. For microparticles, narrow Mie resonances create a strong wavelength dependence of the optical force, enabling more selective sorting. Several beam parameters, such as power, wavelength, polarization state and waveguide modes can be tuned to optimize the sorting performance. The analysis focuses on cylindrical waveguides, where meter-long liquid waveguides in the form of hollow-core photonic crystal fibers are readily available. The modes of such fibers are well-approximated by the cylindrical waveguide modes considered in the theory.

1 Introduction

Since Ashkin’s pioneering paper [1], Optical Tweezers (OT) have opened up the possibility to perform micro- and nanoscale non–contact mechanical manipulation and measurements on particles in the nano-micron scale. Single and multiple optical beams can hold, stretch, and compress particles, filaments, vesicles in a suspension, enabling novel applications such as optical force spectroscopy of microspheres [2, 3]. The room temperature force sensitivity of OT based techniques ranges from 20 fN [4] to 200 pN [5], making this tool ideal to handle particles in the biological micro–environment. Single–beam optical tweezers have allowed the development of several applications in life science, and use optical forces created by a tightly focused laser beam not only to manipulate micron to nanometer–sized structures, but also to measure their mechanical properties [6–8]. The single–beam optical tweezers require gradient forces that exceed scattering forces, a condition usually met with in the case of small differences between medium/particle refractive indexes. This situation is commonly found in aqueous solutions, but not in air or vacuum. The Brownian motion sets the limit for the size of the smallest particle that can be trapped, although at low temperatures (not compatible with live cells), it is possible to even optically trap atomic gases [9]. The success of the single–beam optical tweezers has encouraged the development of new, complementary optical trapping configurations, including optical chromatography, a method that combines the use of optical forces to manipulate and size–select particles in microfluidic systems [6].

Optical chromatography was initially proposed by Imasaka et al. in 1995 [10]. Their approach consisted of mildly focusing a laser beam into a solution containing particles that are counter–flowing coaxially with a capillary. When in equilibrium, there is a balance between the laser scattering force and the drag force from the fluid medium. A challenge in this configuration is to keep the laser beam aligned at the center of the of capillary. Since these pioneering experiments, optical chromatography has led to a range of different applications, for example: to visualize immunological reactions by flowing a sample of antigens past optically–trapped beads that are coated with antibodies. Selective binding changes the effective size of the particles an thus their retention distance. This allowed to determine antigen concentrations in the range from nanomolar to nanograms per milliliter [11]. Optical chromatography can also be used to separate same–sized particles of different refractive indices [12], and be combined with microfluidic channels with well–defined velocity flow patterns [13]. This approach was used to measure low concentrations of anthrax bacteria, by accumulating them from a flow in an optical trap [14]. Further improvements of the technique include hydrodynamic focusing, enabling the separation of equal–sized microsphere differing only 0.1 in refractive index [15]. Further optical chromatography and other optical sorting strategies were summarized in a recent review [6].

Optical chromatography requires a good understanding of the forces acting on the analyzed particles. Initial models were based on geometrical–optics models, and identified analogies with separation parameters for conventional chromatography such as: retention distance, selectivity, dynamic range, theoretical plate number, and resolution [16]. Using a focused Gaussian beam, the retention distance can be related to the Rayleigh length, and can be shown to scale as the square root of particle size. The selectivity in chromatography is directly related to the slope of the retention distance vs. particle size, which can be tuned by carefully changing the laser power and the flow velocity. The dynamic range, defined as the ratio between the maximum and minimum particle size that can be analyzed, is ultimately limited by the size-range of particles that can be trapped. The theoretical plate number, another chromatography parameter, determines the separation efficiency in optical chromatography. In the case of optical traps, this is typically limited by Brownian fluctuations and instrumental noise. Finally, the resolution of particle separation can be determined from the theoretical plate number and selectivity. In theory, particles can be resolved if they differ from at least 0.4% in size [16].

A more recent approach to optical separation uses hollow–core photonic crystal fibers with core diameters between 10 and 30 μm, to optically trap and guide particles along a length of fiber [17]. Advantages are the presence of well–defined waveguide modes, that are well-approximated by the modes of cylindrical waveguides [18], as well as the long interaction length. Optical forces provided by the fundamental fiber mode can robustly align particles in the center of the channel, resulting in well–defined optical- and fluidic forces that can be used to determine the refractive index and size of spherical microspheres by monitoring their speed by Doppler–velocimetry [19–21]. Besides the long distance, stability, and robustness of such waveguide-based chromatography, it also allows to operate in either solutions, air, or vacuum, and at a wide range of pressures and temperatures, as long as the integrity of the waveguide is not compromised. Moreover, the waveguide enclosure around the particles protects the operators and guarantees a perturbation–free environment. Considering the nature of the particles, protection can be an issue in cases such as aerosols loaded with viruses, or radioactive and toxic particles.

The further development of fiber–based optical chromatography techniques requires a better understanding of the optical forces exerted on the particles by the waveguide modes, whose description strongly depends on particles size. The dimensionless size parameter used in scattering is defined as $x=ka=2\pi n_{1}a/\lambda$, where a is the particle’s radius, λ is the wavelength and M = n₁/n₀ is the relative refraction index between the particle and the medium. In the case of Rayleigh particles with radius $x\ll 1$, radial forces are provided by the gradient force which is proportional to x³. If a particle is accelerated in a fiber loop, its maximum velocity will depend on the centripetal force required to keep the particle in its circular trajectory. As both the gradient force and centripetal force scale with x³, the maximum particle velocity will be independent of particle size. The axial forces, on the other hand, are provided by the scattering force and therefore proportional to x⁶. This results in a particle acceleration that is proportional to x³. In the laminar flow regime, the terminal velocity off Rayleigh particles will depend on the ratio between the optical force (x⁶) and the Stokes drag force (x), and will therefore be proportional to x⁵, offering excellent conditions for size–sorting the particles through their velocity. For microparticles, on the other hand, particle resonances and their coupling to the waveguide modes must be taken into account to obtain expressions for the size–dependent optical forces. This requires an integrated scattering theory that combines a Mie description of the particles with analytical waveguide modes, which has not been demonstrated to date.

In this paper, we will develop analytical equations for optical forces, including the scattering forces, on particles within waveguides and will discuss their implications to sorting of particles based on particle size, refraction index, absorption, and mass density. We will start from the full electromagnetic force theory to obtain expressions for the forces acting on waveguide–embedded particles spanning the Rayleigh and Mie regime. Analytical results provide more than just fast calculations compared to numerical methods. They provide powerful insights of how to analyze the role of each parameter, allowing to optimize the design of the system for a given task. The main challenge of analytical optical force calculations for spherical particles lies in the expansion of the fields into vector spherical wave functions, which, except for the plane wave case developed by Mie in 1908 [22], could only be obtained numerically. However, in 2010 we discovered how to perform this expansion for any type of electromagnetic field in the Fourier space and provided analytical integrals for the Beam Shape Coefficients for general waveguides. Analytical results from the integrals were found in the cases of cylindrical and rectangular metallic waveguides [23, 24].

Our current analysis focuses on cylindrical metallic waveguides, whose modes are largely similar to those in readily available hollow-core fibers waveguides. However, the work can be easily extended to rectangular waveguides, using recently develop analytical expressions for their modes [24], which can be of interest for typical on-chip microfluidic geometries. Recently, Ref. [25] presented an analytical solution for the scattered intensities of a sphere on the axis of an infinite metallic cylinder with a plane wave incident on one end of the cylinder. While the combination of the spherical and cylindrical symmetries is interesting, they did not calculate the BSCs for off-axis positions, nor the optical forces on the sphere.

In our analysis, the particle size is assumed to be much smaller than the waveguide size, resulting in a scattering intensity much smaller than the incident light. This implies that we can neglect reflections from the waveguide walls back to the particle, as well as interference between different particles. Furthermore, the waveguides are assumed to be lossless. For low–loss waveguides, losses to be taken into account by considering a slow variation in the electromagnetic field intensity with axial position.

Our manuscript is organized as follows: Section 2 presents the electromagnetic fields for any propagation mode of the cylindrical waveguide. Section 3 defines the Mie Coefficients as well is the Beam Shape Coefficients (BSCs), which are obtained for any waveguide mode and any particle position. These results are valid for any particle size smaller than the waveguide size and any waveguide mode. In Section 5 the BSCs will be used to obtain a complete theory for the optical forces, for any waveguide mode, and any particle size and position within the waveguide. At this point the theory is complete.

We will then discuss the implication of the theory for three specific applications. First we focus on nanoparticles that are well in the Rayleigh regime. We will obtain their optical forces using a small–value series expansion of Mie Coefficients, and will demonstrate that radial forces are proportional to x³, and axial forces proportional to x⁶ for non–absorbing particles. Secondly, we discuss the microparticles, for which the full Mie calculation is required. Section 6, studies the optical force as a function of particle size and refractive index in case of an on-axis particle position. For larger particles, strong Mie resonances appear that can potentially be used be used to perform in highly–sensitive particle sorting experiments. Third, in Section 7, we will study how the effect of particle absorption on axial force depends on particle size. Finally, we will summarize illustrative results for special cases, applications, and experimental implementation.

2 Waveguides Electromagnetic Fields

The general formalism for the TM and TE modes of a hollow waveguide of arbitrary cross–section starts with the scalar solution $g\left({\mathbf{r}}^{\prime }\right)=g_{\gamma }\left({\mathit{\rho }}^{\prime }\right)e^{\text{i}{k}_z{z}^{\prime }}$ of the Helmholtz equation in three dimensions

$\left[\frac{{\text{d}}^2}{\text{d}{\mathbf{r}}^{\prime }2}+k^2\right]g\left({\mathbf{r}}^{\prime }\right)=0,(1)$

which leads to a two–dimensional solution $g_{\gamma }\left(\mathit{\rho }\right)$ of the transverse Helmholtz equation [26].

$\left[\frac{{\text{d}}^2}{\text{d}{\mathit{\rho }}^{\prime }2}+{\gamma }^2\right]g_{\gamma }\left({\mathit{\rho }}^{\prime }\right)=0,(2)$

satisfying the boundary conditions at the waveguide surface S. The TM modes satisfy $g_{\gamma }\left(\mathit{\rho }\right){|}_S=0$ and the TE modes satisfy ${\partial }_n{g}_{\gamma }\left(\mathit{\rho }\right){|}_S=0$ with ${\partial }_n$ the normal derivative, k_z is the wavevector in the z–direction, $k^2=k_z^2+{\gamma }^2$, γ is the transverse wavevector, $\text{d}/\text{d}\mathit{\rho }=\widehat{x}{\partial }_x+\widehat{y}{\partial }_y$ is the transverse gradient operator, and $\text{d}/\text{d}\mathbf{r}=\text{d}/\text{d}\mathit{\rho }+\widehat{z}{\partial }_z$ is the three–dimensional gradient operator, $k^2=n_o{\omega }^2/c^2$, n_o is the waveguide’s interior refractive index. It is useful to introduce the cylindrical coordinate system in the r–space and k–space. In the r–space we have $\mathit{\rho }=x\widehat{x}+y\widehat{y}$, ${\rho }^2=\mathit{\rho }\cdot \mathit{\rho }$ and $\widehat{\mathit{\rho }}=\mathit{\rho }/\rho$. We have also $\varphi =-y\widehat{x}+x\widehat{y}$ and $\widehat{\varphi }=\varphi /\rho$. In the same way, for spherical coordinates in k–space, we change $\left(r,\theta ,\varphi \right)\to \left(k,\xi ,\zeta \right)$ and for cylindrical coordinates, we change $\left(\rho ,\varphi ,z\right)\to \left(\gamma ,\zeta ,k_z\right)$, obtaining equivalent systems of coordinates. The boundary conditions define possible values for γ based on the geometry alone. The cutoff frequency of a given mode happens when $k^2={\gamma }^2$. Below this frequency, k_z < 0, the mode becomes evanescent. Given the scalar solution, the TM and TE fields are given by:

$\left[\begin{array}{c}{\mathbf{E}}^{\text{TM}}\\ Z{\mathbf{H}}^{\text{TE}}\end{array}\right]=\frac{k_z}{k}\widehat{z}\times \left[\begin{array}{c}-Z{\mathbf{H}}^{\text{TM}}\\ {\mathbf{E}}^{\text{TE}}\end{array}\right]=E_0[\widehat{z}+\text{i}\frac{k_z}{{\gamma }^2}\frac{\text{d}}{\text{d}{\mathit{\rho }}^{\prime }}]g\left({\mathbf{r}}^{\prime }\right).(3)$

For a general waveguide, the task is to find g(r) satisfying boundary conditions. For a cylindrical waveguide, the scalar wavefunction is given by $g_{\gamma }\left(\mathit{\rho }\right)=J_m\left({\gamma }_{mn}\rho \right)e^{\pm \text{i}m\varphi }$, where $J_m(.)$ are Bessel functions. The transverse vectors are given by ${\gamma }_{mn}^{\text{TE}}={\chi }_{mn}^{\prime }$ where ${\chi }_{mn}^{\prime }$ is the nth–root of $J_m(.)$ derivative, that is $J_m^{\prime }\left({\chi }_{mn}^{\prime }\right)=0$ for the TE mode, and ${\gamma }_{mn}^{\text{TM}}={\chi }_{mn}$ where ${\chi }_{mn}$ is the n^th–root of $J_m(.)$, that is $J_m\left({\chi }_{mn}\right)=0$, for the TM mode. The cutoff frequencies are ${\omega }_{cmn}^{\text{TE}}=\left(c/n_o\right)\left({\chi }_{mn}^{\prime }/R\right)$ for TE modes, and ${\omega }_{cmn}^{\text{TM}}=\left(c/n_o\right)\left({\chi }_{nm}/R\right)$ for TM mode. The cutoff wavelengths are given by ${\lambda }_{cnm}=2\pi R/{\chi }_{nm}$. The first ten modes, from highest to lowest ${\lambda }_{nm}$ are ${\text{TE}}_{11}\left[{\chi }_{11}^{\prime }\right]=1.8412$, ${\text{TM}}_{10}\left[{\chi }_{10}\right]=2.4048$, ${\text{TE}}_{21}\left[{\chi }_{21}^{\prime }\right]=3.0542$, ${\text{TE}}_{01}\left[{\chi }_{01}^{\prime }\right]=3.8317$, ${\text{TM}}_{11}\left[{\chi }_{11}\right]=3.8317$, ${\text{TE}}_{31}\left[{\chi }_{31}^{\prime }\right]=4.2012$, ${\text{TM}}_{12}\left[{\chi }_{12}\right]=5.1356$, ${\text{TE}}_{41}\left[{\chi }_{41}^{\prime }\right]=5.3175$, ${\text{TE}}_{12}\left[{\chi }_{21}^{\prime }\right]=5.3314$, ${\text{TM}}_{20}\left[{\chi }_{20}\right]=5.5201$. Even for $R=2.5\mu \text{m}$ and n_o = 1 there are 16 modes for wavelengths shorter than 2,000 nm and 45 modes for wavelengths shorter than 1,000 nm. Therefore, in the visible–NIR region, the waveguide will be essentially multimode. We are mostly interested in the range 400 to 2,000 nm to avoid direct IR and NIR absorption. Assuming a wavelength of 500 nm and a $10\mu \text{m}$ waveguide diameter, the first mode’s propagation constant becomes $k_z/k=0.9995\approx 1$ and $\gamma /k=\sqrt{1-{\left(k_z/k\right)}^2}=0.032$. Linear combinations of $g^{\pm }\left(\mathit{\rho }\right)$ can be used to build the cosine and sine modes, as well.

The analytical expression for the fields of TE, TM, ±, cos and sin are condensed in the expressions:

$\frac{{\mathbf{E}}^{\text{TM}}}{E_0{e}^{\text{i}{k}_{z}z}}=\left(\text{i}\frac{k_z}{\gamma }{J}_m^{\prime }\left(\gamma \rho \right)\widehat{\mathit{\rho }}+J_m\left(\gamma \rho \right)\widehat{z}\right)\left[\begin{array}{c}{e}^{\pm \text{i}m\varphi }\\ \cos m\varphi \\ \sin m\varphi \end{array}\right]-\frac{k_z}{\gamma }\frac{m}{\gamma \rho }{J}_m\left(\gamma \rho \right)\widehat{\varphi }\left[\begin{array}{c}\pm e^{\pm \text{i}m\varphi }\\ \text{i}\sin m\varphi \\ -\text{i}\cos m\varphi \end{array}\right](4)$

$\frac{Z{\mathbf{H}}^{\text{TM}}}{E_0{e}^{\text{i}{k}_{z}z}}=\frac{k}{\gamma }\frac{m}{\gamma \rho }{J}_m\left(\gamma \rho \right)\widehat{\mathit{\rho }}\left[\begin{array}{c}\pm e^{\pm \text{i}m\varphi }\\ \text{i}\sin m\varphi \\ -\text{i}\cos m\varphi \end{array}\right]+\text{i}\frac{k}{\gamma }{J}_m^{\prime }\left(\gamma \rho \right)\widehat{\varphi }\left[\begin{array}{c}{e}^{\pm \text{i}m\varphi }\\ \cos m\varphi \\ \sin m\varphi \end{array}\right].(5)$

where the form TE/TM and 3 × 1 array represents the (±, cos, sin) forms and $Z=\sqrt{\mu /\mathit{ϵ}}$ is the waveguide medium impedance. It is interesting to note the interchange between electric and magnetic fields between TM and TE modes does not change the result. This shows that the Poynting vector expression will be the same for both the TE and TM modes. The TE mode is obtained by: ${\mathbf{E}}^{\text{TE}}\to -Z{\mathbf{H}}^{\text{TM}}$ and $Z{\mathbf{H}}^{\text{TE}}\to {\mathbf{E}}^{\text{TM}}$. The first modes have m = 1 and for a very large waveguide ($R\to \infty$, tending to free space), $\gamma \to 0$ and $k_z\to k$, we can use $\gamma \rho \to 0$ to obtain $J_m\left(x\right)/x\approx x^m/\left(2^{m}m!\right)$, $J_m^{\prime }\left(x\right)\approx m{x}^{m-1}/2^{m}m!$ and so $J_1\left(x\right)/x=J_1^{\prime }\left(x\right)=1/2$, which allows us to obtain:

$\frac{{\mathbf{E}}^{\text{TM}}}{E_0}=\frac{i}{2}\frac{k_z}{k}\left[\begin{array}{c}\widehat{x}\pm \text{i}\widehat{y}\\ \widehat{x}\\ \widehat{y}\end{array}\right]e^{ikz}.(6)$

This is clearly the plane wave limit, from which we can notice the association of circularly polarized light for the mode ±, x-polarized beam to cos mode and y-polarized to sin mode. One can couple to ± modes using circular polarized light while a polarized beam will couple to cos and sin modes.

We also have the explicit form of the electromagnetic fields written as

$|\frac{{\mathbf{E}}^{\text{TE}}}{E_0}{|}^2=\frac{k^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2{J}_m^2\left(\gamma \rho \right)\left[\begin{array}{c}1\\ {\text{sin}}^{2}m\varphi \\ {\text{cos}}^{2}m\varphi \end{array}\right]+\frac{k^2}{{\gamma }^2}{J}_m^{\prime }|\left(\gamma \rho \right)\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right](7)$

$|\frac{Z{\mathbf{E}}^{\text{TM}}}{E_0}{|}^2=\frac{k_z^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2{J}_m^2\left(\gamma \rho \right)\left[\begin{array}{c}1\\ {\text{sin}}^{2}m\varphi \\ {\text{cos}}^{2}m\varphi \end{array}\right]+(\frac{k_z^2}{{\gamma }^2}{J}_m^{\prime }\left(\gamma \rho \right)+J_m^2\left(\gamma \rho \right))\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right],(8)$

and it is important to calculate the gradients, given by

$\begin{array}{c}{\nabla }_{\rho }|\frac{{\mathbf{E}}^{\text{TE}}}{E_0}{|}^2=2\gamma \frac{k^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2\left(J_m\left(\gamma \rho \right)J_m^{\prime }\left(\gamma \rho \right)-\frac{1}{\gamma \rho }{J}_m^2\left(\gamma \rho \right)\right)\left[\begin{array}{c}1\\ {\text{sin}}^{2}m\varphi \\ {\text{cos}}^{2}m\varphi \end{array}\right]+\\ 2\gamma \frac{k^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2{J}_m^{\prime }\left(\gamma \rho \right)J_m^{″}\left(\gamma \rho \right)\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right]\end{array}(9)$

and

$\begin{array}{c}{\nabla }_{\rho }|\frac{{\mathbf{E}}^{\text{TM}}}{E_0}{|}^2=2\gamma \frac{k^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2\left(J_m\left(\gamma \rho \right)J_m^{\prime }\left(\gamma \rho \right)-\frac{1}{\gamma \rho }{J}_m^2\left(\gamma \rho \right)\right)\left[\begin{array}{c}1\\ {\text{sin}}^{2}m\varphi \\ {\text{cos}}^{2}m\varphi \end{array}\right]+\\ 2\gamma \left(\frac{k^2}{{\gamma }^2}{\left(\frac{m}{\gamma \rho }\right)}^2{J}_m^{\prime }\left(\gamma \rho \right)J_m^{″}\left(\gamma \rho \right)+J_m\left(\gamma \rho \right)J_m^{\prime }\left(\gamma \rho \right)\right)\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right].\end{array}(10)$

The ϕ component of the gradient is given by

${\nabla }_{\varphi }|\frac{{\mathbf{E}}^{\text{TM}}}{E_0}{|}^2=\gamma \frac{k^2}{{\gamma }^2}\frac{1}{\gamma \rho }({\left(\frac{m}{\gamma \rho }{J}_m\left(\gamma \rho \right)\right)}^2-J_m^{\prime }2\left(\gamma \rho \right))\left[\begin{array}{c}0\\ 1\\ -1\end{array}\right](11)$

and

$\begin{array}{c}{\nabla }_{\varphi }|\frac{{\mathbf{E}}^{\text{TM}}}{E_0}{|}^2=\frac{\gamma }{\gamma \rho }\frac{k_z^2}{{\gamma }^2}({\left(\frac{m}{\gamma \rho }{J}_m\left(\gamma \rho \right)\right)}^2-J_m^{\prime }|\left(\gamma \rho \right))\left[\begin{array}{c}0\\ -1\\ 1\end{array}\right]+\\ \frac{\gamma }{\gamma \rho }{J}_m^2\left(\gamma \rho \right)\left[\begin{array}{c}0\\ 1\\ -1\end{array}\right]\end{array}(12)$

3 Vector Spherical Wave Function Expansion

To calculate the electromagnetic fields inside and scattered by a spherical particle located at an arbitrary position r inside a cylindrical waveguide, one has to expand the incident fields in the Hansen multipoles. For the force calculation, only the incident and scattered fields are necessary. The VSWF expansion is given by

$\left[\begin{array}{c}{\mathbf{E}}^{\text{in}}\left({\mathbf{r}}^{\prime }\right)\\ Z{\mathbf{H}}^{\text{in}}\left({\mathbf{r}}^{\prime }\right)\end{array}\right]=E_0\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\left(\left[\begin{array}{c}{G}_{pq}^{\text{TE}}\\ G_{pq}^{\text{TM}}\end{array}\right]{\mathbf{M}}_{pq}^{\text{in}}\left({\mathbf{r}}^{\prime }\right)+\left[\begin{array}{c}{G}_{pq}^{\text{TM}}\\ -G_{pq}^{\text{TE}}\end{array}\right]{\mathbf{N}}_{pq}^{\text{in}}\left({\mathbf{r}}^{\prime }\right)\right)(13)$

$\left[\begin{array}{c}{\mathbf{E}}^{\text{sc}}\left({\mathbf{r}}^{\prime }\right)\\ Z{\mathbf{H}}^{\text{sc}}\left({\mathbf{r}}^{\prime }\right)\end{array}\right]=-E_0\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\left(\left[\begin{array}{c}{b}_p{G}_{pq}^{\text{TE}}\\ a_p{G}_{pq}^{\text{TM}}\end{array}\right]{\mathbf{M}}_{pq}^{\text{sc}}\left({\mathbf{r}}^{\prime }\right)+\left[\begin{array}{c}{a}_p{G}_{pq}^{\text{TM}}\\ -b_p{G}_{pq}^{\text{TE}}\end{array}\right]{\mathbf{N}}_{pq}^{\text{sc}}\left({\mathbf{r}}^{\prime }\right)\right),(14)$

where E₀ is an electric field dimension constant, $k{\mathbf{N}}_{pq}^{\text{in}}\left(\mathbf{r}\right)=\text{i}\nabla \times {\mathbf{M}}_{pq}^{\text{in}}\left(\mathbf{r}\right)$, and ${\mathbf{M}}_{pq}^{\text{in}}\left(\mathbf{r}\right)=j_p\left(kr\right){\mathbf{X}}_{pq}\left(\widehat{r}\right)$ are the Hansen multipoles for the incident beam, and $k{\mathbf{N}}_{pq}^{\text{sc}}\left(\mathbf{r}\right)=\text{i}\nabla \times {\mathbf{M}}_{pq}^{\text{sc}}\left(\mathbf{r}\right)$, and ${\mathbf{M}}_{pq}^{\text{sc}}\left(\mathbf{r}\right)=h_p^{\left(1\right)}\left(kr\right){\mathbf{X}}_{pq}\left(\widehat{r}\right)$ the Hansen multipoles for the scattered beam, j_p(kr) are spherical Bessel functions, y_p(kr) are spherical Neumann functions and $h_p^{\left(1\right)}\left(kr\right)=j_p\left(kr\right)+\text{i}{y}_p\left(kr\right)$${\mathbf{X}}_{pq}\left(\widehat{r}\right)=\mathbf{L}{Y}_{pq}\left(\widehat{r}\right)/\sqrt{p\left(p+1\right)}$ are the vector spherical harmonics, $Y_{pq}\left(\widehat{r}\right)$ are the scalar spherical harmonics, $Z=\sqrt{\mu /\epsilon }$, $k=\omega \sqrt{\mu \epsilon }$, and $\mathbf{L}=-i\mathbf{r}\times \text{d}/\text{d}\mathbf{r}$ is the angular momentum operator in direct space.

The coefficients $G_{pq}^{\text{TE}/\text{TM}}$ are known as Beam Shape Coefficients (BSC) [27]. It is important to note that, as the VSWF expansion is performed with the origin of the axis at the center of a sphere, we need to apply boundary conditions that are independent of the particular position of the beam. The a_p and b_p coefficients, called Mie coefficients, come from the boundary conditions over the surface of a sphere, and are given by

$a_p=\frac{M{\psi }_p^{\prime }\left(x\right){\psi }_p\left(Mx\right)-{\psi }_p\left(x\right){\psi }_p^{\prime }\left(Mx\right)}{M{\xi }_p^{\prime }\left(x\right){\psi }_p\left(Mx\right)-{\xi }_p\left(x\right){\psi }_p^{\prime }\left(Mx\right)}(15)$

$b_p=\frac{M{\psi }_p\left(x\right){\psi }_p^{\prime }\left(Mx\right)-{\psi }_p^{\prime }\left(x\right){\psi }_p\left(Mx\right)}{M{\xi }_p\left(x\right){\psi }_p^{\prime }\left(Mx\right)-{\xi }_p^{\prime }\left(x\right){\psi }_p\left(Mx\right)}.(16)$

Here $x=ka=2\pi n_{1}a/\lambda$ is the size parameter, where a is the particle’s radius and λ is the wavelength and M = n₁/n₀ is the relative refraction index. The functions $\psi \left(x\right)=x{j}_p\left(x\right)$ and $\xi \left(x\right)=x{h}_p^{\left(1\right)}\left(x\right)$ are known as Riccati-Bessel functions. The BSCs depend on the incident field only and all further calculations rely on them. The position of the center of the sphere with respect to the axis of the incident field in cylindrical coordinates is defined as $\mathbf{r}=\left({\rho }_o,{\varphi }_o,z_o\right)$. It is clear, then, that BSCs are functions of this position, $G_{pq}^{\text{TE},\text{TM}}=G_{pq}^{\text{TE},\text{TM}}\left(\mathbf{r}\right)$. Gustav Mie solved this problem in 1908 [22] for the incident plane wave using the Rayleigh expansion of a plane wave in terms of spherical Bessel functions. For the plane wave, the position of the sphere with respect to the beam is not important. For an arbitrary beam, the usual procedure is to take the dot product with ${\mathbf{X}}_{pq}\left(\widehat{r}\right)$ to obtain:

$E_0{j}_p\left(k{r}^{\prime }\right)\left[\begin{array}{c}{G}_{pq}^{\text{TE}}\\ G_{pq}^{\text{TM}}\end{array}\right]={{\displaystyle \int }}^{\text{}}\text{dΩ}\left(\widehat{r}\prime \right){\mathbf{X}}_{pq}^{\text{*}}\left(\widehat{r}\prime \right)\cdot \left[\begin{array}{c}\mathbf{E}\left({\mathbf{r}}^{\prime }\right)\\ Z\mathbf{H}\left({\mathbf{r}}^{\prime }\right)\end{array}\right],(17)$

which is not complete because there is a spherical Bessel function embedded in the right–hand side that did not cancel explicitly the one on the left side. We solved this problem for an arbitrary electromagnetic field in 2010 [23, 24], by taking the Fourier transform, explicitly canceling the spherical Bessel functions, and obtaining the final result:

$\left[\begin{array}{c}{G}_{pq}^{\text{TE}}\\ G_{pq}^{\text{TM}}\end{array}\right]=\frac{{\text{i}}^p}{E_0}\sqrt{\frac{2}{\pi }}{{\displaystyle \int }}^{\text{}}\text{dΩ}\left(\widehat{k}\prime \right){\mathrm{}}_{pq}^{\text{*}}\left(\widehat{k}\prime \right)\cdot \left[\begin{array}{c}{\mathrm{ℰ}}_{\text{k}}\left({\widehat{k}}^{\prime }\right)\\ Z{\mathrm{\mathscr{H}}}_{\text{k}}\left({\widehat{k}}^{\prime }\right)\end{array}\right],(18)$

where

$\left[\begin{array}{c}\mathbf{E}\left({\mathbf{r}}^{\prime }\right)\\ \mathbf{H}\left({\mathbf{r}}^{\prime }\right)\end{array}\right]=\frac{1}{{\left(2\pi \right)}^{3/2}}{{\displaystyle \int }}^{\text{}}{\text{d}}^3{k}^{\prime }\left[\begin{array}{c}\mathrm{ℰ}\left({\mathbf{k}}^{\prime }\right)\\ \mathrm{\mathscr{H}}\left({\mathbf{k}}^{\prime }\right)\end{array}\right]e^{\text{i}\mathbf{k}\prime \cdot \mathbf{r}}(19)$

are Fourier transforms of each other and

$\left[\begin{array}{c}\mathrm{ℰ}(\mathbf{k}\prime )\\ \mathrm{\mathscr{H}}(\mathbf{k}\prime )\end{array}\right]=\frac{\delta \left(k^{\prime }-k\right)}{k^{\prime 2}}\left[\begin{array}{c}{\mathrm{ℰ}}_{\text{k}}\left({\widehat{k}}^{\prime }\right)\\ {\mathrm{\mathscr{H}}}_{\text{k}}\left({\widehat{k}}^{\prime }\right)\end{array}\right].(20)$

are the Fourier transforms in a sphere, requiring that $|\mathbf{k}\prime |=k$. By using the translation theorem of Fourier transform it was easy to calculate the BSC for any position respect to the incident beam. The geometry of the problem is shown in Figure 1.

FIGURE 1

FIGURE 1. Cylindrical waveguide coordinate system (CWG). The natural coordinate system of the waveguide is rʺ and a point r at which we perform a partial wave expansion (PWE) in vector spherical wave functions (VSWF). In the new system, the expansion is done in the rʹ coordinate system, the fields are described by the Hansen multipoles and the G^TE/TM coefficients, that are constant in the rʹ frame, but varies depending on the position of the expansion r.

4 Cylindrical Waveguide BSC’s

Using Eqs 24, 25 and 32 of Ref. [24] the explicit expression for the BSCs for any mode of the waveguide is given by:

$\begin{array}{c}\left[\begin{array}{cc}{G}_{pq}^{\text{TE}}\left[\text{TE}\right]& G_{pq}^{\text{TE}}\left[\text{TM}\right]\\ G_{pq}^{\text{TM}}\left[\text{TE}\right]& G_{pq}^{\text{TM}}\left[\text{TM}\right]\end{array}\right]={\text{i}}^{p-q}\sqrt{\frac{4\pi \left(2p+1\right)}{p\left(p+1\right)}\frac{\left(p-q\right)!}{\left(p+q\right)!}}\left[\begin{array}{cc}\text{i}{P}_p^q\prime \left(\frac{k_z}{k}\right)& {\left(\frac{k}{\gamma }\right)}^{2}q{P}_p^q\left(\frac{k_z}{k}\right)\\ {\left(\frac{k}{\gamma }\right)}^{2}q{P}_p^q\left(\frac{k_z}{k}\right)& -\text{i}{P}_p^q\prime \left(\frac{k_z}{k}\right)\end{array}\right]\\ \times e^{-\text{i}q{\varphi }_o}{e}^{\text{i}{k}_z{z}_o}{g}_{mq}^{\text{form}}({\rho }_o).\end{array}(21)$

Here, our notation is that superscript TE/TM means the spherical Transverse Electric/Magnetic in the sphere, while (TE)/(TM) means the waveguide TE/TM modes, $P_p^q(.)$ is the associated Legendre function. To shortening the equations we will omit the arguments with the convention that all Legendre functions $P_p^q=P_p^q\left(k_z/k\right)$, and all Bessel functions $J_m=J_m\left(\gamma {\rho }_o\right)$. The form factor is:

$g_{m,q}^{\text{form}}\left[\begin{array}{c}\pm \\ \text{cos}\\ \text{sin}\end{array}\right]=\left[\begin{array}{c}{(\pm )}^q{J}_{m\mp q}{e}^{\pm \text{i}m{\varphi }_o}\\ \frac{1}{2}\left(J_{m-q}+{\left(-1\right)}^q{J}_{m+q}\right)\text{cos}m{\varphi }_o+\text{i}\frac{1}{2}\left(J_{m-q}-{\left(-1\right)}^q{J}_{m+q}\right)\text{sin}m{\varphi }_o\\ \frac{1}{2}\left(J_{m-q}+{\left(-1\right)}^q{J}_{m+q}\right)\text{sin}m{\varphi }_o-\text{i}\frac{1}{2}\left(J_{m-q}-{\left(-1\right)}^q{J}_{m+q}\right)\text{cos}m{\varphi }_o\end{array}\right].(22)$

5 Optical Forces

In a recent review, it was demonstrated how optical forces can be derived from the Maxwell stress tensor given the BSCs in the Minkowski formalism [28]. Explicitly, in Eq. 93 in Ref. [28] provides the expression for the transverse forces, while the force in the z–direction is given by Eq. 94 in Ref. [28]:

$\begin{array}{c}{F}_z=-\frac{\epsilon |E_o{|}^2}{2k^2} \mathfrak{Re}{\displaystyle \sum }_{pq}\frac{\text{i}}{p+1}\{\sqrt{\frac{p\left(p+2\right)\left(p+q+1\right)\left(p-q+1\right)}{\left(2p+1\right)\left(2p+3\right)}}\\ \times [A_p{G}_{p+1,q}^{\text{TM}}{G}_{p,q}^{\text{TM*}}+B_p{G}_{p+1,q}^{\text{TE}}{G}_{p,q}^{\text{TE*}}]-\frac{q}{p}{C}_p{G}_{p,q}^{\text{TM}}{G}_{p,q}^{\text{TE*}}\}\end{array}(23)$

where $A_p=a_p+a_{p+1}^{\text{∗}}-2a_p^{\text{∗}}{a}_{p+1}$, $B_p=b_p+b_{p+1}^{\text{∗}}-2b_p^{\text{∗}}{b}_{p+1}$ and $C_p=a_p+b_p^{\text{∗}}-2a_p^{\text{∗}}{a}_p^{\text{∗}}$. The definition of the literature, including the a recent review [28], is $A_p^{lit}=A_p^{\text{*}}$ which is counter-intuitive, because the lowest order Mie coefficient appears conjugated $a_p^{\text{*}}$. We therefore prefer the current definition. The formulas must thus be corrected by changing $A_p^{lit}\to A_p^{\text{*}}$ and $B_p^{lit}\to B_p^{\text{*}}$ and vice-versa, while keeping the same C_p.

$\begin{array}{c}\left[\begin{array}{c}{F}_{\rho }\\ F_{\varphi }\end{array}\right]=-\frac{\epsilon |E_o{|}^2}{{\left(2k\right)}^2}\left[\begin{array}{c}\mathfrak{Re}\\ \mathfrak{Im}\end{array}\right]{\displaystyle \sum }_{p=1}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\frac{\text{i}{e}^{-\text{i}\varphi }}{p+1}\{\sqrt{\frac{p\left(p+2\right)\left(p+q+2\right)\left(p+q+1\right)}{\left(2p+1\right)\left(2p+3\right)}}\\ \times [A_p{G}_{p+1,-\left(q+1\right)}^{\text{TM}}{G}_{p,-q}^{\text{TM*}}+A_p^{\text{*}}{G}_{p,q}^{\text{TM}}{G}_{p+1,q+1}^{\text{TM*}}+B_p{G}_{p+1,-\left(q+1\right)}^{\text{TE}}{G}_{p,-q}^{\text{TE*}}+B_p^{\text{*}}{G}_{p,q}^{\text{TE}}{G}_{p+1,q+1}^{\text{TE*}}]+\\ -\frac{1}{p}\sqrt{\left(p-q\right)\left(p+q+1\right)}[C_p{G}_{p,q}^{\text{TM}}{G}_{p,q+1}^{\text{TE*}}-C_p^{\text{*}}{G}_{p,q}^{\text{TE}}{G}_{p,q+1}^{\text{TM*}}]\}\end{array}(24)$

From the expression of z-component of the force we substitute the BSCs expressions into the force equation, and using the normalization for the force $\mathbf{C}=\mathbf{F}/\left(\pi \epsilon |E_o/k{|}^2\right)$ we obtain

$\begin{array}{l}\frac{1}{2}{C}_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\mathfrak{Re}\left\{A_p^{\text{*}}\right\}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\frac{\left(p+q+1\right)}{{\left(p+1\right)}^2}\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}\frac{k^4}{{\gamma }^4}{q}^2{P}_p^q{P}_{p+1}^q\\ P_p^{\prime q}{P}_{p+1}^q\prime \end{array}\right]|g_{m,q}^{\text{form}}{|}^2+\\ \stackrel{\infty }{{\displaystyle \sum }_{p=1}}\mathfrak{Re}\left\{B_p^{\text{*}}\right\}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\frac{\left(p+q+1\right)}{{\left(p+1\right)}^2}\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}{P}_p^{\prime q}{P}_{p+1}^q\prime \\ \frac{k^4}{{\gamma }^4}{q}^2{P}_p^q{P}_{p+1}^q\end{array}\right]|g_{m,q}^{\text{form}}{|}^2\\ -\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\mathfrak{Re}\left\{C_p^{\text{*}}\right\}\stackrel{p}{{\displaystyle \sum }_{q=-p}}\frac{\left(2p+1\right)}{p^2{\left(p+1\right)}^2}\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}1\\ 1\end{array}\right]\frac{k^2}{{\gamma }^2}{q}^2{P}_p^q{P}_p^q\prime |g_{m,q}^{\text{form}}{|}^2.\end{array}(25)$

We display the equation in this way to demonstrate that all particle properties, such as size and refractive index, are embedded in the Mie factors A_p, B_p and C_p. Breaking the q sum into q = 0, q > 0 and q < 0, and using $\left(p+q\right)!P_p^{-q}={\left(-1\right)}^q\left(p-q\right)!P_p^q$, we can simplify this expression to:

$\begin{array}{l}\frac{1}{2}{C}_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\frac{\mathfrak{Re}\left\{A_p^{\text{*}}\right\}}{\left(p+1\right)}\left\{\left[\begin{array}{c}0\\ P_p{\prime }^0{P}_{p+1}^0\prime \end{array}\right]\right|g_{m,0}^{\text{form}}{|}^2\\ +\stackrel{p}{{\displaystyle \sum }_{q=1}}\frac{1}{\left(p+1\right)}\frac{\left(p-q+1\right)!}{\left(p+q\right)!}\left[\begin{array}{c}\frac{k^4}{{\gamma }^4}{q}^2{P}_p^q{P}_{p+1}^q\\ P_p{\prime }^q{P}_{p+1}^q\prime \end{array}\right]\left(\right|g_{m,q}^{\text{form}}{|}^2+\left|g_{m,-q}^{\text{form}}{|}^2\right)\}\\ +\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\frac{\mathfrak{Re}\left\{B_p^{\text{*}}\right\}}{\left(p+1\right)}\left\{\left[\begin{array}{c}{P}_p{\prime }^0{P}_{p+1}^0\prime \\ 0\end{array}\right]\right|g_{m,0}^{\text{form}}{|}^2\\ +\stackrel{p}{{\displaystyle \sum }_{q=1}}\frac{1}{\left(p+1\right)}\frac{\left(p-q+1\right)!}{\left(p+q\right)!}\left[\begin{array}{c}{P}_p{\prime }^q{P}_{p+1}^q\prime \\ \frac{k^4}{{\gamma }^4}{q}^2{P}_p^q{P}_{p+1}^q\end{array}\right]\left(\right|g_{m,q}^{\text{form}}{|}^2+\left|g_{m,-q}^{\text{form}}{|}^2\right)\}\\ -\stackrel{\infty }{{\displaystyle \sum }_{p=1}} \mathfrak{Re}\left\{C_p^{\text{*}}\right\}\stackrel{p}{{\displaystyle \sum }_{q=1}}\frac{\left(2p+1\right)}{p^2{\left(p+1\right)}^2}\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}1\\ 1\end{array}\right]\frac{k^2}{{\gamma }^2}{q}^2{P}_p^q{P}_p^q\prime \\ \times \left(\right|g_{m,q}^{\text{form}}{|}^2+\left|g_{m,-q}^{\text{form}}{|}^2\right).\end{array}(26)$

The [radial, azimuthal] forces can be written in the same way:

$\begin{array}{c}\frac{1}{2}{C}_{[\rho ,\varphi ]}\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]={\displaystyle \sum _{p=1}^{\infty }\frac{\mathfrak{I}\mathfrak{m}\left\{A_p\right\}}{{\left(p+1\right)}^2}\left\{{\displaystyle \sum _{q=1}^p\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}\frac{k^4}{{\gamma }^4}q\left(q+1\right)P_p^q{P}_{p+1}^{q+1}\\ {P^{\prime }}_p^q{P}_{p+1}^{q+1}\prime \end{array}\right]\times \left[\mathfrak{R}\mathfrak{e},\mathfrak{I}\mathfrak{m}\right]\left\{g_{m,q}^{\text{form}}{g}_{m,q+1}^{form\text{*}}-g_{m,q}^{\text{form}}{g}_{m,-q-1}^{form\text{*}}\right\}}\right\}}\\ +{\displaystyle \sum _{p=1}^{\infty }\frac{\mathfrak{I}\mathfrak{m}\left\{B_p\right\}}{{\left(p+1\right)}^2}\left\{{\displaystyle \sum _{q=1}^p\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}{P^{\prime }}_p^q{P}_{p+1}^{q+1}\prime \\ \frac{k^4}{{\gamma }^4}q\left(q+1\right)P_p^q{P}_{p+1}^{q+1}\end{array}\right]\times \left[\mathfrak{R}\mathfrak{e},\mathfrak{I}\mathfrak{m}\right]\left\{g_{m,q}^{\text{form}}{g}_{m,q+1}^{form\text{*}}-g_{m,q}^{\text{form}}{g}_{m,-q-1}^{form\text{*}}\right\}}\right\}}\\ +{\displaystyle \sum _{p=1}^{\infty }\frac{\mathfrak{R}\mathfrak{e}\left\{A_p\right\}}{{\left(p+1\right)}^2}\left\{{\displaystyle \sum _{q=1}^p\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}\frac{k^4}{{\gamma }^4}q\left(q+1\right)P_p^q{P}_{p+1}^{q+1}\\ {P^{\prime }}_p^q{P}_{p+1}^{q+1}\prime \end{array}\right]\times \left[\mathfrak{I}\mathfrak{m},\mathfrak{R}\mathfrak{e}\right]\left\{g_{m,q}^{\text{form}}{g}_{m,q+1}^{form\text{*}}-g_{m,q}^{\text{form}}{g}_{m,-q-1}^{form\text{*}}\right\}}\right\}}\\ +{\displaystyle \sum _{p=1}^{\infty }\frac{\mathfrak{R}\mathfrak{e}\left\{B_p\right\}}{{\left(p+1\right)}^2}\left\{{\displaystyle \sum _{q=1}^p\frac{\left(p-q\right)!}{\left(p+q\right)!}\left[\begin{array}{c}{P^{\prime }}_p^q{P}_{p+1}^{q+1}\prime \\ \frac{k^4}{{\gamma }^4}q\left(q+1\right)P_p^q{P}_{p+1}^{q+1}\end{array}\right]\times \left[\mathfrak{I}\mathfrak{m},\mathfrak{R}\mathfrak{e}\right]\left\{g_{m,q}^{\text{form}}{g}_{m,q+1}^{form\text{*}}-g_{m,q}^{\text{form}}{g}_{m,-q-1}^{form\text{*}}\right\}}\right\}}\\ -{\displaystyle \sum _{p=1}^{\infty }\frac{\mathfrak{R}\mathfrak{e}\left\{C_p\right\}}{{\left(p+1\right)}^2}\left\{{\displaystyle \sum _{q=1}^p\frac{2p+1}{p^2}\frac{\left(p-q\right)!}{\left(p+q\right)!}\frac{k^2}{{\gamma }^2}\left[q{P}_p^q{P^{\prime }}_p^{q+1}-\left(q+1\right){P^{\prime }}_p^q{P}_p^{q+1}\right]\times \left[\mathfrak{I}\mathfrak{m},\mathfrak{R}\mathfrak{e}\right]\left\{g_{m,q}^{\text{form}}{g}_{m,q+1}^{form\text{*}}\right\}\left[\begin{array}{c}1\\ -1\end{array}\right]}\right\}}\end{array}(27)$

While these equations could be written in a more compact way, we want to highlight the dependence of the real and imaginary parts of the Mie factors A_p, B_p and C_p, because these factors define the Mie resonances, as well as their scaling with particle size and relative refraction index, which are important to the sorting process. The dependence of the position of the particle inside the waveguide is embedded in the form factors $g_{m,q}^{\text{form}}\left({\rho }_o,{\varphi }_o\right)$.

It becomes immediately clear that F_z depends on the real part of A_p, B_p and C_p only. In principle, $F_{\rho }$ depends on the real and imaginary parts of A_p, B_p and C_p, but the form factors will cancel the real part. It is also interesting to see that dependence on z disappeared, as expected from the axial translation symmetry of the waveguide (assumed lossless). Furthermore, there is no azimuthal dependence for the ± modes, corresponding to circular polarization, due to their rotational symmetry. For the $\text{cos}\varphi$ and $\text{sin}\varphi$ forms, the result does depend on ϕ, as the symmetry is broken by the incident beam polarization.

These expressions are general and can be used to obtain the force for any waveguide mode and particle size. Special care must be taken in the choice of highest p to truncate the series. The Wiscombe criterion was to truncate the sum when $|a_p{|}^2+|b_p{|}^2<5{10}^{-14}$ [29], leading to the known truncation of the infinite series of $p_{\text{max}}=\text{int}\left\{x+4.05x^{1/2}+2\right\}$. Meanwhile, Neves and Pisignano included the relative error, allowing faster computation within an fixed error bound, as $p_{\text{max}}=\text{int}\left\{x+0.76{\mathit{ϵ}}^{2/3}{x}^{1/3}-4.1\right\}$, where ϵ is related to a relative error of ${10}^{-\mathit{ϵ}}$ [30].

In the following sections we will discuss two special cases. The first is related to use of the optical chromatography (OC) for nanoparticles that are well in the Rayleigh limit $x\ll 1$. In this limit, we need the Mie coefficients at $x\to 0$. Our results will show that for any waveguide mode in any position inside the guide, the radial forces are proportional to x³ and the gradient of the field intensity $F_{\rho ,\varphi }\propto x^3\nabla |\mathbf{E}{|}^2$. On the other hand, for non–absorbing particle, the axial force is proportional to x⁶ and the field intensity $F_z\propto x^6|\mathbf{E}{|}^2$. This shows that radial forces are gradient forces while z forces are scattering forces. The axial force result is obvious because there is no intensity gradient in a lossless waveguide.

The second case is waveguide-based OC for microparticles, such as dielectric spheres or virus–loaded droplets, with sizes in the Mie regime $x\approx 1$. In the Mie–limit, the radial force is assumed to keep the particle on the z-axis, $\rho =0$, resulting in a significant simplification of the expressions. The analysis can be done for any particle size and any waveguide mode, allowing us to observe the role played by the Mie resonances in the optical sorting process.

6 Mie Coefficients in Rayleigh Regime

We can use the series expansion of spherical Bessel and Hankel function keeping the lowest orders of the imaginary and real parts to obtain [31]:

$a_p\approx -\text{i}\frac{\left(p+1\right)\left(M^2-1\right)x^{2p+1}}{\left(2p+1\right)!!\left(2p-1\right)!!\left(p+1+p{M}^2\right)}+{\left[\frac{\left(p+1\right)\left(M^2-1\right)x^{2p+1}}{\left(2p+1\right)!!\left(2p-1\right)!!\left(p+1+p{M}^2\right)}\right]}^2(28)$

$b_p\approx -\text{i}\frac{\left(M^2-1\right)x^{2p+3}}{\left(2p+3\right)!!\left(2p+1\right)!!}+{\left[\frac{\left(M^2-1\right)x^{2p+3}}{\left(2p+3\right)!!\left(2p+1\right)!!}\right]}^2(29)$

For the lowest terms we obtain

$a_1\approx -\text{i}\left[\frac{2}{3}\frac{M^2-1}{M^2+2}\right]x^3+{\left[\frac{2}{3}\frac{M^2-1}{M^2+2}\right]}^2{x}^6(30)$

$b_1\approx -\text{i}\left[\frac{M^2-1}{45}\right]x^5+{\left[\frac{M^2-1}{45}\right]}^2{x}^{10}(31)$

Reports on the scattering of metallic nanoparticles usually call the lowest terms approximations as modified long–wavelength approximation [32–35]. We shall not neglect the second term because they are real while the first are imaginary. Except for TE[m = 0] modes, we will keep only terms with $a_1$ meaning that transverse forces scale with x³ while the axial force scales with x⁶. However, for particles with absorption, the refractive index becomes complex in the form $M=M^{\prime }+\text{i}{M}^{\prime \prime }$. In this case the lowest order terms are:

$a_1\approx -\text{i}\mathfrak{Re}\left\{\frac{2}{3}\frac{M^2-1}{M^2+2}\right\}{x}^3+\mathfrak{Im}\left\{\frac{2}{3}\frac{M^2-1}{M^2+2}\right\}{x}^3(32)$

$b_1\approx -\text{i}\mathfrak{Re}\left\{\frac{M^2-1}{45}\right\}{x}^5+\mathfrak{Im}\left\{\frac{M^2-1}{45}\right\}{x}^5(33)$

Without absorption, the leading term comes only from a₁, and only A₁ = a₁ and C₁ = a₁ shall be considered. The only $P_p^q$ to take into account are given by $P_2^2=3{\gamma }^2/k^2$, $P_2^1=-3\gamma k_z/k^2$ and $P_2^0=1/2\left[\left(3k_z^2/k^2\right)-1\right]$ for the q =2 family whose derivative are $P_2^2\prime =-6\left(kz/k\right)$, $P_2^1\prime =3\left(k_z^2/k\gamma -\gamma /k\right)$, $P_2^0\prime =3\left(k_z/k\right)$, and for q = 1 $P_1^1=-\gamma /k$, $P_1^1\prime =k_z/\gamma$, $P_1^0=k_z/k$, $P_1^0\prime =1$.

In the expressions for the force, we get

$C_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=3\frac{k_z}{k}\mathfrak{Re}\left\{a_1\right\}\left(2\frac{{\gamma }^2}{k^2}|g_{m,0}^{\text{form}}{|}^2\left[\begin{array}{c}0\\ 1\end{array}\right]+\left(|g_{m,1}^{\text{form}}{|}^2+|g_{m,-1}^{\text{form}}{|}^2\right)\left[\begin{array}{c}1\\ k_z^2/k^2\end{array}\right]\right)(34)$

for the z-component and

$\begin{array}{c}{C}_{\rho }\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\frac{3}{2}\frac{k}{\gamma } \mathfrak{Im}\{-a_1\}( \mathfrak{Re}\{g_{m,1}^{\text{form}}{g}_{m,2}^{\text{form}}-g_{m,-2}^{\text{form}}{g}_{m,-1}^{\text{form}}\}\left[\begin{array}{c}1\\ k_z^2/k^2\end{array}\right]+\\ + \mathfrak{Re}\{g_{m,0}^{\text{form}}{g}_{m,1}^{\text{form}}-g_{m,-1}^{\text{form}}{g}_{m,0}^{\text{form}}\}\left[\begin{array}{c}-1\\ 2{\gamma }^2/k^2+k_z^2/k^2\end{array}\right])\\ -\frac{3}{2}\frac{k}{\gamma } \mathfrak{Re}\left\{a_1\right\}(\mathfrak{Im}\{g_{m,1}^{\text{form}}{g}_{m,2}^{\text{form}}+g_{m,-2}^{\text{form}}{g}_{m,-1}^{\text{form}}\}\left[\begin{array}{c}1\\ k_z^2/k^2\end{array}\right]\\ + \mathfrak{Re}\{g_{m,0}^{\text{form}}{g}_{m,1}^{\text{form}}+g_{m,-1}^{\text{form}}{g}_{m,0}^{\text{form}}\}\left[\begin{array}{c}-1\\ 2{\gamma }^2/k^2+k_z^2/k^2\end{array}\right])\end{array}(35)$

for the radial component, while the azimuthal one is given by

$\begin{array}{c}{C}_{\varphi }\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\frac{3}{2}\frac{k}{\gamma } \mathfrak{Im}\{-a_1\}(\mathfrak{Im}\{g_{m,1}^{\text{form}}{g}_{m,2}^{\text{form}}-g_{m,-2}^{\text{form}}{g}_{m,-1}^{\text{form}}\}\left[\begin{array}{c}1\\ k_z^2/k^2\end{array}\right]+\\ + \mathfrak{Im}\{g_{m,0}^{\text{form}}{g}_{m,1}^{\text{form}}-g_{m,-1}^{\text{form}}{g}_{m,0}^{\text{form}}\}\left[\begin{array}{c}-1\\ 2{\gamma }^2/k^2+k_z^2/k^2\end{array}\right])\\ +\frac{3}{2}\frac{k}{\gamma } \mathfrak{Re}\left\{a_1\right\}(\mathfrak{Re}\{g_{m,1}^{\text{form}}{g}_{m,2}^{\text{form}}+g_{m,-2}^{\text{form}}{g}_{m,-1}^{\text{form}}\}\left[\begin{array}{c}1\\ k_z^2/k^2\end{array}\right]+\\ + \mathfrak{Re}\{g_{m,0}^{\text{form}}{g}_{m,1}^{\text{form}}+g_{m,-1}^{\text{form}}{g}_{m,0}^{\text{form}}\}\left[\begin{array}{c}-1\\ 2{\gamma }^2/k^2+k_z^2/k^2\end{array}\right])\end{array}(36)$

Next the form factors are calculated. We can bring all of them in terms of J_m(x) instead of $J_{m\mp 1}\left(x\right)$ and $J_{m\mp 2}\left(x\right)$ combinations. Using the Bessel function differential equation $J_m^{\prime \prime }\left(x\right)=-J_m^{\prime }/x-\left(1-m^2/x^2\right)J_m$ and the canonical three term recurrence relations, we can calculate the above combinations of $g_{m,q}^{\text{form}}$. The first thing one can see is that $\mathfrak{Im}\left\{g_{m,0}^{\text{form}}{g}_{m,1}^{\text{form}}+g_{m,-1}^{\text{form}}{g}_{m,0}^{\text{form}}\right\}=0$ and $\mathfrak{Im}\left\{g_{m,1}^{\text{form}}{g}_{m,2}^{\text{form}}+g_{m,2}^{\text{form}}{g}_{m,1}^{\text{form}}\right\}=0$, because these sums are always real. Using $u=\gamma {\rho }_o$ to shortening of the equations in algebraic manipulations, one can promptly find that

$F_z^{TE}=\frac{3\pi \epsilon }{k^2}\frac{k_z}{k}\mathfrak{Re}\left\{a_1\right\}|{\mathbf{E}}^{TE}{|}^2(37)$

$=\frac{3\pi \epsilon }{{\gamma }^2}\frac{k_z}{k}|E_o{|}^2\mathfrak{Re}\left\{a_1\right\}\left[\begin{array}{c}{J}_m^{\prime }\left(u\right)+{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right)\\ J_m^{\prime }\left(u\right){\text{cos}}^{2}m\varphi +{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right){\text{sin}}^{2}m\varphi \\ J_m^{\prime }\left(u\right){\text{sin}}^{2}m\varphi +{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right){\text{cos}}^{2}m\varphi \end{array}\right](38)$

in the three line one column notation for the form factors. The TM case is then written as

$F_z^{TM}=\frac{3\pi \epsilon }{k^2}\frac{k_z}{k}\mathfrak{Re}\left\{a_1\right\}|{\mathbf{E}}^{TM}{|}^2(39)$

$=\frac{3\pi \epsilon }{{\gamma }^2}\frac{k_z}{k}|E_o{|}^2\mathfrak{Im}\left\{a_1\right\}\frac{k_z^2}{k^2}\left[\begin{array}{c}{J}_m^2\left(u\right)+{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right)\\ J_m^2\left(u\right){\text{cos}}^{2}m\varphi +{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right){\text{sin}}^{2}m\varphi \\ J_m^2\left(u\right){\text{sin}}^{2}m\varphi +{\left(\frac{m}{u}\right)}^2{J}_m^2\left(u\right){\text{cos}}^{2}m\varphi \end{array}\right]+\frac{3\pi \epsilon }{{\gamma }^2}\frac{k_z}{k}|E_o{|}^2\mathfrak{Im}\left\{a_1\right\}\frac{{\gamma }^2}{k^2}{J}_m^2\left(u\right)\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right](40)$

This means that

$F_z=\frac{4\pi \epsilon k_z}{k^3}|{\mathbf{E}}^{TM}{|}^2{\left(\frac{M^2-1}{M^2+2}\right)}^2{x}^6(41)$

which is proportional to particle radius to the sixth power and to the wavelength to the inverse forth power, as expected int he Rayleigh regime. The radial force in Rayleigh regime is given by

$F_{\rho }^{TE}=\frac{3\pi \epsilon }{2k^3}\mathfrak{Im}\left\{-a_1\right\}\nabla \rho |{\mathbf{E}}^{TM}{|}^2(42)$

$C_{\rho }^{TE}=6\frac{k}{\gamma }\mathfrak{Im}\left\{-a_1\right\}\left[\begin{array}{c}{J}_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right)+\frac{m^2}{u^2}\left(J_m\left(u\right)J_m^{\prime }\left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right)\\ J_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right){\text{cos}}^{2}m\varphi +\frac{m^2}{u^2}\left(J_m\left(u\right)J_m^{\prime }\left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right){\text{sin}}^{2}m\varphi \\ J_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right){\text{sin}}^{2}m\varphi +\frac{m^2}{u^2}\left(J_m\left(u\right)J_m^{\prime }\left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right){\text{cos}}^{2}m\varphi \end{array}\right],(43)$

and

$F_{\rho }^{TM}=\frac{3\pi \epsilon }{2k^3}\mathfrak{Im}\left\{-a_1\right\}\nabla \rho |{\mathbf{E}}^{TM}{|}^2(44)$

$C_{\rho }^{TM}=6\frac{k_z}{\gamma }\mathfrak{Im}\left\{-a_1\right\}\left[\begin{array}{c}{J}_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right)+\frac{m^2}{u^2}\left(J_m\left(u\right)J_m\prime \left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right)\\ J_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right){\text{cos}}^{2}m\varphi +\frac{m^2}{u^2}\left(J_m\left(u\right)J_m\prime \left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right){\text{sin}}^{2}m\varphi \\ J_m^{\prime \prime }\left(u\right)J_m^{\prime }\left(u\right){\text{sin}}^{2}m\varphi +\frac{m^2}{u^2}\left(J_m\left(u\right)J_m\prime \left(u\right)-\frac{1}{u}{J}_m^2\left(u\right)\right){\text{cos}}^{2}m\varphi \end{array}\right]+6\frac{\gamma }{k}\mathfrak{Re}\left\{-a_1\right\}{J}_m^{\prime }\left(u\right)J_m\left(u\right)\left[\begin{array}{c}1\\ {\text{cos}}^{2}m\varphi \\ {\text{sin}}^{2}m\varphi \end{array}\right].(45)$

Therefore, the radial force

$F_{\rho }=\frac{\pi \epsilon }{k^3}{\nabla }_{\rho }|\mathbf{E}{|}^2\left(\frac{M^2-1}{M^2+2}\right)x^3(46)$

is proportional to the particle’s radius to the third power. The maximum radial force exceeds the axial scattering force, and keeps the particle near the position of maximum of the field intensity. Although some waveguide modes a zero intensity on axis, the commonly excited fundamental waveguide mode does tend to bring particles back to it's center. Finally, we discuss the azimuthal force. This force has two terms, one proportional to $\mathfrak{Im}\left\{-a_1\right\}$, a gradient force, an the other proportional to $\mathfrak{Re}\left\{a_1\right\}$. By taking the lowest order term, with $\mathfrak{Im}\left\{-a_1\right\}$, we obtain

$F_{\varphi }^{TE}=\frac{\pi \epsilon }{k^3}{\nabla }_{\varphi }|\mathbf{E}{|}^2\left(\frac{M^2-1}{M^2+2}\right)x^3(47)$

$C_{\varphi }^{TE}=3\frac{k}{\gamma } \mathfrak{Im}\{-a_1\}\frac{1}{u}(J_m^{\prime }\left(u\right)-\frac{m^2}{u^2}{J}_m^2\left(u\right))\left[\begin{array}{c}-0\\ -1\\ 1\end{array}\right]2m\sin m\varphi \cos m\varphi (48)$

that shows a gradient force type for the transverse force, and scattering force type for the axial force,

$\left[\begin{array}{c}{F}_{\rho }\\ F_{\varphi }\\ F_z\end{array}\right]=\epsilon \frac{M^2-1}{M^2+2}{R}^3\left[\begin{array}{c}\pi {\nabla }_{\rho }|\mathbf{E}{|}^2\\ \pi {\nabla }_{\varphi }|\mathbf{E}{|}^2\\ \frac{M^2-1}{M^2+2}\gamma k^3\left(\frac{4}{3}\pi R^3\right)|\mathbf{E}{|}^2\end{array}\right].(49)$

7 Forces Acting on-Axis Particles

For on-axis Mie-particles, ${\rho }_o=0$, the only Bessel function that survives are the $J_0\left(0\right)=1$ with $J_p\left(0\right)={\delta }_{p,0}$. As expected by symmetry, the transverse forces are zero because there is a mix between $J_{m\pm q}$ with $J_{m\pm q\pm 1}$. Only the axial force is different from zero. The form factors become

$|g_{m,q}^{\text{form}}{|}^2+|g_{m,-q}^{\text{form}}{|}^2={\delta }_{m,q}\left[\begin{array}{c}1\\ 1/2\\ 1/2\end{array}\right],(50)$

which cancels the sum in q, and

$|g_{m,0}^{\text{form}}{|}^2={\delta }_{m,0}\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right].(51)$

Clearly the sin form is null for m = 0 and cos form is just half of ±. The result ± form is

$\begin{array}{c}\frac{1}{4}{C}_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]={\delta }_{m,0}\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\frac{P_p{\prime }^0{P}_{p+1}^0\prime }{\left(p+1\right)}\left[\begin{array}{c}\mathfrak{Re}\left\{B_p\right\}\\ \mathfrak{Re}\left\{A_p\right\}\end{array}\right]\\ +\left(1-{\delta }_{m,0}\right)\stackrel{\infty }{{\displaystyle \sum }_{p=m}}\frac{\mathfrak{Re}\left\{A_p\right\}}{{\left(p+1\right)}^2}\frac{\left(p-m+1\right)!}{\left(p+m\right)!}\left[\begin{array}{c}\frac{k^4}{{\gamma }^4}{m}^2{P}_p^m{P}_{p+1}^m\\ P_p^{\prime m}{P}_{p+1}^m\prime \end{array}\right]\\ +\left(1-{\delta }_{m,0}\right)\stackrel{\infty }{{\displaystyle \sum }_{p=m}}\frac{\mathfrak{Re}\left\{B_p\right\}}{{\left(p+1\right)}^2}\frac{\left(p-m+1\right)!}{\left(p+m\right)!}\left[\begin{array}{c}{P}_p^{\prime m}{P}_{p+1}^m\prime \\ \frac{k^4}{{\gamma }^4}{m}^2{P}_p^m{P}_{p+1}^m\end{array}\right]\\ -\left(1-{\delta }_{m,0}\right)\stackrel{\infty }{{\displaystyle \sum }_{p=m}}\frac{\mathfrak{Re}\left\{C_p\right\}}{{\left(p+1\right)}^2}\frac{\left(2p+1\right)}{p^2}\frac{\left(p-m\right)!}{\left(p+m\right)!}\left[\begin{array}{c}1\\ 1\end{array}\right]\frac{k^2}{{\gamma }^2}{m}^2{P}_p^m{P}_p^m\prime \end{array}(52)$

For the first modes where $k_z/k\approx 1$ and $\gamma /k=\sqrt{1-{\left(k_z/k\right)}^2}$ we can use the fact that $\sqrt{1-x^2}{P}_p^m\prime \left(x\to 1\right)={\delta }_{m,1}p\left(p+1\right)/2$, $P_p^0\prime \left(x\to 1\right)=p\left(p+1\right)/2$ and $P_p^m\left(1\right)={\delta }_{m,0}$ to simplify the force to

$\frac{1}{4}{C}_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\stackrel{\infty }{{\displaystyle \sum }_{p=1}}p\left(p+2\right)\left(\left(p+1\right){\delta }_{m,0}+\frac{k^2}{{\gamma }^2}\frac{{\delta }_{m,1}}{p+1}\right)\left[\begin{array}{c}\mathfrak{Re}\left\{A_p\right\}\\ \mathfrak{Re}\left\{B_p\right\}\end{array}\right](53)$

which are nonzero only for the m = 0 and m = 1 modes. Also, for m = 0 modes the force for any particle size simplifies to:

$\frac{1}{4}{C}_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\stackrel{\infty }{{\displaystyle \sum }_{p=1}}\frac{P_p^0\prime P_{p+1}^0\prime }{p+1}\left[\begin{array}{c}\mathfrak{Re}\left\{A_p\right\}\\ \mathfrak{Re}\left\{B_p\right\}\end{array}\right].(54)$

On axis force in the Rayleigh regime can be obtained by making p = 1, $\mathfrak{Re}\left\{A_p\right\}= \mathfrak{Re}\left\{C_p\right\}= \mathfrak{Re}\left\{a_1\right\}$ and $\mathfrak{Re}\{B_p=0\}$. Now, $0\le m\le 1$, due to the Legendre Functions, therefore only modes with m = 0 and m = 1 exists, and the force is given by

$C_z\left[\begin{array}{c}\text{TE}\\ \text{TM}\end{array}\right]=\frac{4}{3}\frac{k_z}{k}\left(\frac{M^2-1}{M^2+2}\right)x^6\left(2{\delta }_{m,0}\left[\begin{array}{c}0\\ 1\end{array}\right]+{\delta }_{m,1}\left[\begin{array}{c}{k}^2/{\gamma }^2\\ k_z^2/{\gamma }^2\end{array}\right]\right)(55)$

8 Results

We analyzed the important role of Mie Resonances, for different refractive index particles as a function of their sizes. Typical glass particles and liquid droplets have refractive indices in the range between 1.3 and 1.5, but semiconductors can have n as high as 3. Usually the high refractive index appears at the absorption band edge, except for quantum dots (QD), for which the optical gap is larger than the bulk band gap. Although QD are typically smaller than 10 nm sizes, one could also consider highly QD–doped microparticles with optical properties that are dominated by the QD refractive index. Therefore, in this analysis, we considered two limits of 1.3 and 3.0 for the refractive index. To exemplify this, the axial force for TE₁₀ and TM₁₀ modes was calculated as a function of particle radius for a n = 1.3 and n = 3.0 refractive index, assuming a air-filled waveguide with a diameter of $R=10\mu m$, see Figure 2.

FIGURE 2

FIGURE 2. Axial force for: left n = 1.3 particle; right n = 3 particle. Blue TE₁₀ mode, red TM₁₀ mode of a 10 μm radius cylindrical waveguide.

For both particles, the resonances becoming narrower as the radius increases, with the narrowest resonances occurring for the n = 3 particle. The optical force landscape of the two particle’s parameters, size and relative refractive index is shown in map plot in Figure 3. The results suggest that Mie resonances can indeed be used to discriminate between a narrow range of discrete particles.

FIGURE 3

FIGURE 3. Density plot of axial force as function of particle’s relative refractive index M and particle radius for a 10 μm radius cylindrical waveguide.

Figures 4 and 5 show the forces for the first four first TE and TM modes for a Rayleigh particle inside a cylindrical waveguide. The axial force F_z is clearly proportional to the local field intensity $|\mathbf{E}{|}^2$, and $F_{\rho }$ to ${\nabla }_{\rho }|\mathbf{E}{|}^2$. The modes that allows keep the particle on the axis are that ones that have energy in concentrated in the axis such as the TE₁₁ mode or the TM₁₁ mode, for which $F_{\rho }<0$, providing a returning force. This relation depends on the value of $k/\gamma$ as can be seen in Eqs. 7 and 8. It can also be noted that, in the Rayleigh regime, $F_{\rho }\gg F_z$.

FIGURE 4

FIGURE 4. First TE modes: optical forces seen as profile from the center to the walls of the waveguide. Depending on the mode $F_{\rho }$ can be always positive, always negative, or change the sign across the waveguide profile. Calculations were performed using the full expressions on Eq 5, assuming a waveguide radius of 5 µm, and a size parameter x = 0.1. Note that $F_z\propto |\mathbf{E}{|}^2$.

FIGURE 5

FIGURE 5. First TM modes: optical forces seen as profile from the center to the walls of the waveguide. As we have the $\widehat{\mathit{\rho }}$ pointing outwards the waveguide profile, the $F_{\rho }<0$ acts as contrary to this direction, bringing the particle to the center of the waveguide. Calculations were performed using the full expressions on Eqs 23 and 24, assuming a waveguide radius of 5 µm, and a size parameter x = 0.1.

In the above calculation, we have used the full expressions for the force calculation. We used a wavelength of $\lambda =1064\text{ nm}$, a waveguide radius of 5 µm, and index of refraction of the internal medium of 1.01 and a particle’s index of refraction of 1.33 which gives as relative index of refraction M = 1.317. The particle radius was chosen to be 16.76 nm, resulting in a size parameter of x = 0.1. This gives a truncation criterium of l_MAX = 4, according to [30]. The values of the wavevector parameters are $k=5.96\times {10}^6 {\text{m}}^{-1}$, $k_z=5.94\times {10}^6 {\text{m}}^{-1}$, $\gamma =4.81\times {10}^5{\text{ m}}^{-1}$, and so $k_z/k=0.9967$ and $\gamma /k=0.08$. The other values are the roots of Bessel functions and of their derivatives, and universal constants $\epsilon =8.854\times {10}^{-12}\text{ F}/\text{m}$, $\mu =1.256\times {10}^{-6}\text{ N}/{\text{A}}^2$ both in SI units, and Z = 376.639 SI units.

Figure 6 shows how the force varies with size parameter (and in this case is directly proportional to the size of the particle x = ka). It also shows the truncation criteria for calculating these forces. In this case, the forces are of the same magnitude (we are far from the Rayleigh regime). There are some challenges to use this methodology to perform particle sorting. The first is how to load the particles into the waveguide, to start the sorting process.

FIGURE 6

FIGURE 6. At the top the coefficients for x = 20.8751 (for a particle of 3.5 µm and a refractive index of 3. At bottom the forces (from left right F_z, $F_{\rho }$ and $F_{\varphi }$) as function of the size parameter (on the legend at right), from the lower one (particle of 1.5 µm) to the biggest one (particle of radius 3.5 µm). We changed the radius of the waveguide to 20 µm to allow for larger values of x.

9 Conclusion and Discussion

We have demonstrated analytical expressions for optical forces acting on particles inside cylindrical waveguides. The theory is valid for particle radius a ranging from the Rayleigh regime to large microparticles, and for arbitrary particle position and optical modes. Analytical results obtained for a range of different particles sizes and refractive indices demonstrate that Mie-resonances can result in a strong size-dependence of the scattering force, with implications for waveguide-based particle sorting. In the Rayleigh regime, the radial, axial, and azimuthal forces were studied as a function of particle position within the waveguide. In the Rayleigh regime, the most important factor for particles sorting is found to be the x³, x⁶ size scale of the optical forces, followed by the refractive index differences, and the presence of absorption.

A possible application of this work would be post–synthesis size selection of quantum dots. Quantum dots are nanoparticles in the 1–10 nm diameter sizes, well in the Rayleigh regime, which show a size dependent optical band gap. We have produced quantum dots by chemical synthesis and, also, by Pulsed Laser Deposition (PLD) [36–38]. The advantage of PLD is a total degree of freedom regarding the material, just by focusing a powerful laser beam on any piece of solid bulk material, to produce a broad range of nano/micro-particles. It can be used with radioactive materials, semiconductors, or oxides, which can be hard to produce by chemical synthesis. PLD, however, produces a much wider size distribution than the chemical synthesis methods. Hence, the possibility to sort the quantum dots can make this technique more attractive, especially inside a waveguide where no particles are lost and that can be reprocessed. For this application, one can use the size–dependent absorption edges. With a tunable laser one can start with the longer wavelengths for which only the largest QDs absorb. Because the force in the presence of absorption scales with x³, much larger than the x⁶ non-absorbing QDs, one can isolate biggest absorbing QDs. The wavelength of the laser can then be swept, from the long to short, to perform an absorption–selective size separation. QDs usually have cap layers that change their hydrodynamic radii and their diffusion times in a liquid medium [39]. It would therefore also be possible to QDs of equal semiconductor core radius, through their cap layer thickness.

Furthermore, for microparticles ($2\pi x>\lambda$), particles exhibit strong Mie resonances, [9, 10] that can be used to perform a fine selection of particle sizes. Experiments on virus-loaded erosols have become an important topic in the current covid-19 pandemic, because the air human-human disease transmission is defined mainly by the sizes of the droplets. Large droplets fall to the ground fast, while smaller ones can fly much farther away. Because the hollow-core fiber waveguides have a transparent glass cladding, one could even measure the amount of virus in each particle using a “optical flow cytometry” approach.

For solid particles, we suggest to follow the particle–loading process used by [20], which mechanically removes the particles from a glass plate, and uses two counter-propagation beams to confine them to a region in front of the waveguide entrance. To launch a particles, the counter propagating beam is switched off and the particle is optically guided through the waveguide. The collection and sorting at the end of the waveguide can be done by applying a sweeping transverse force to the particles the end of the waveguide, that directs the articles to a specific region on the collector - analogous to a mass spectrometer. There are a range of parameters that can be adjusted for optimal sorting performance, including wavelength, polarization, and optical power, all of which can be changed fast. In case of charges particles, commonly found in colloids and cells, transverse electric fields could be applied to the waveguide [40] to selectively collect or slow down particles based on their charge.

Experimentally available hollow waveguides are typically based on hollow-core photonic crystal fibers, for which the theory described here offers a reasonable approximation but no analytical expression. It will therefore be interesting to explore different ways to create small-diameter cylindrical hollow metallic waveguides. One approach could be based on coating metal on the interior of a tube of a material with a thermal softening point close to that of the metal and subsequently draw this to optical fiber. Before the drawing, some propagation modes could be suppressed by laser cutting the metal perpendicular to the wall’s currents of the targeted modes.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

CC, TE, and AN started the development of the theory, never published before, back in 2008 in Max Planck Institute in Erlangen, Germany, and conceived the idea for this report. WM developed the Fourier Transform method to calculate the BSCs. AF and AN performed the calculations of Mie resonances on optical forces. All the authors wrote, edited, reviewed and revised the manuscript. CC supervised the study.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

AF acknowledges the National Institute of Photonics (INFo). CC acknowledges National Institute of Science and Technology in Photonics Applied to Cell Biology [INFABIC] CNPq (573913/2008-00/FAPESP (08/57906-3); CNPq grant (312049/2014-5); “Física do Petróleo em Meios Porosos”, PETROBRAS-UFC (2016/00328-4) and F020/WIPPS II—Simulação numérica de invasão de água em poços produtores”, PETROGAL-UFC T.G.E. acknowledges the support from the Winton Program for the Physics of Sustainability, the Isaac Newton Trust, and the Leverhulme Trust (RPG‐2018‐256).

References

1. Ashkin A. Acceleration and trapping of particles by radiation pressure. Phys Rev Lett (1970) 24:156–9. doi:10.1103/PhysRevLett.24.156

CrossRef Full Text | Google Scholar

2. Fontes A, Neves AAR, Moreira WL, de Thomaz AA, Barbosa LC, Cesar CL, et al. Double optical tweezers for ultrasensitive force spectroscopy in microsphere mie scattering. Appl Phys Lett (2005) 87:221109. doi:10.1063/1.2137896

CrossRef Full Text | Google Scholar

3. Neves AAR, Fontes A, de Y Pozzo L, de Thomaz AA, Chillce E, Rodriguez E, et al. Electromagnetic forces for an arbitrary optical trapping of a spherical dielectric. Optic Express (2006) 14:13101–6. doi:10.1364/OE.14.013101

CrossRef Full Text | Google Scholar

4. Rohrbach A. Switching and measuring a force of 25 femtonewtons with an optical trap. Optic Express (2005) 13:9695–701. doi:10.1364/OPEX.13.009695

CrossRef Full Text | Google Scholar

5. Callegari A, Mijalkov M, Gököz AB, Volpe G. Computational toolbox for optical tweezers in geometrical optics. J Opt Soc Am B (2015) 32:B11–9. doi:10.1364/JOSAB.32.000B11

CrossRef Full Text | Google Scholar

6. Jonáš A, Zemánek P. Light at work: the use of optical forces for particle manipulation, sorting, and analysis. Electrophoresis (2008) 29:4813–51. doi:10.1002/elps.200800484 |

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Fontes A, Fernandes HP, de Thomaz AA, Barbosa LC, Barjas-Castro ML, Cesar CL. Measuring electrical and mechanical properties of red blood cells with double optical tweezers. J Biomed Optic (2008) 13:1–6. doi:10.1117/1.2870108

CrossRef Full Text | Google Scholar

8. Fontes A, Ajito K, Neves AAR, Moreira WL, de Thomaz AA, Barbosa LC, et al. Raman, hyper-Raman, hyper-Rayleigh, two-photon luminescence and morphology-dependent resonance modes in a single optical tweezers system. Phys Rev (2005b) 72:012903. doi:10.1103/PhysRevE.72.012903 |

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Chu S, Bjorkholm JE, Ashkin A, Cable A. Experimental observation of optically trapped atoms. Phys Rev Lett (1986) 57:314–7. doi:10.1103/PhysRevLett.57.314 |

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Imasaka T, Kawabata Y, Kaneta T, Ishidzu Y. Optical chromatography. Anal Chem (1995) 67:1763–5. doi:10.1021/ac00107a003

CrossRef Full Text | Google Scholar

11. Hatano T, Kaneta T, Imasaka T. Application of optical chromatography to immunoassay. Anal Chem (1997) 69:2711–5. doi:10.1021/ac970081q |

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hart SJ, Terray AV. Refractive-index-driven separation of colloidal polymer particles using optical chromatography. Appl Phys Lett (2003) 83:5316–8. doi:10.1063/1.1635984

CrossRef Full Text | Google Scholar

13. Terray A, Arnold J, Hart SJ. Enhanced optical chromatography in a pdms microfluidic system. Optic Express (2005) 13:10406–15. doi:10.1364/OPEX.13.010406

CrossRef Full Text | Google Scholar

14. Hart SJ, Terray A, Arnold J, Leski TA. Sample concentration using optical chromatography. Optic Express (2007) 15:2724–31. doi:10.1364/OE.15.002724 |

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Terray A, Hebert CG, Hart SJ. Optical chromatographic sample separation of hydrodynamically focused mixtures. Biomicrofluidics (2014) 8:064102. doi:10.1063/1.4901824 |

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kaneta T, Ishidzu Y, Mishima N, Imasaka T. Theory of optical chromatography. Anal Chem (1997) 69:2701–10. doi:10.1021/ac970079z |

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Euser TG, Garbos MK, Chen JSY, Russell PS. Precise balancing of viscous and radiation forces on a particle in liquid-filled photonic bandgap fiber. Opt Lett (2009) 34:3674–6. doi:10.1364/OL.34.003674 |

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Marcatili EAJ, Schmeltzer RA. Hollow metallic and dielectric waveguides for long distance optical transmission and lasers. Bell Syst Tech J (1964) 43:1783–809. doi:10.1002/j.1538-7305.1964.tb04108.x

CrossRef Full Text | Google Scholar

19. Garbos MK, Euser TG, Schmidt OA, Unterkofler S, Russell PS. Doppler velocimetry on microparticles trapped and propelled by laser light in liquid-filled photonic crystal fiber. Opt Lett (2011) 36:2020–2. doi:10.1364/OL.36.002020 |

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Schmidt OA, Garbos MK, Euser TG, Russell PSJ. Metrology of laser-guided particles in air-filled hollow-core photonic crystal fiber. Opt Lett (2012) 37:91–3. doi:10.1364/OL.37.000091 |

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Sharma A, Xie S, Zeltner R, Russell PS. On-the-fly particle metrology in hollow-core photonic crystal fibre. Optic Express (2019) 27:34496–504. doi:10.1364/OE.27.034496

CrossRef Full Text | Google Scholar

22. Mie G. Beiträge zur optik trüber medien, speziell kolloidaler metallösungen. Ann Phys (1908) 330:377–445.

CrossRef Full Text | Google Scholar

23. Moreira WL, Neves AAR, Garbos MK, Euser TG, Russell PSJ, Cesar CL. Expansion of arbitrary electromagnetic fields in terms of vector spherical wave functions (2010). e-printsArXiv: 1003.2392.

CrossRef Full Text | Google Scholar

24. Moreira WL, Neves AAR, Garbos MK, Euser TG, Cesar CL. Expansion of arbitrary electromagnetic fields in terms of vector spherical wave functions. Optic Express (2016) 24:2370–82. doi:10.1364/OE.24.002370

CrossRef Full Text | Google Scholar

25. Mangini F, Tedeschi N. Scattering of an electromagnetic plane wave by a sphere embedded in a cylinder. J Opt Soc Am A (2017) 34:760–9. doi:10.1364/JOSAA.34.000760

CrossRef Full Text | Google Scholar

26. Jackson JD. Classical electrodynamics. New York, NY: Wiley (1999).

CrossRef Full Text | Google Scholar

27. Gouesbet G, Gréhan G. Generalized Lorenz-Mie theories. Berlin, Germany: Springer-Verlag (2011).

CrossRef Full Text | Google Scholar

28. Neves AAR, Cesar CL. Analytical calculation of optical forces on spherical particles in optical tweezers: tutorial. J Opt Soc Am B (2019) 36:1525–37. doi:10.1364/JOSAB.36.001525

CrossRef Full Text | Google Scholar

29. Wiscombe WJ. Improved mie scattering algorithms. Appl Optic (1980) 19:1505–9. doi:10.1364/AO.19.001505

CrossRef Full Text | Google Scholar

30. Neves AAR, Pisignano D. Effect of finite terms on the truncation error of mie series. Opt Lett (2012) 37:2418–20. doi:10.1364/OL.37.002418 |

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bohren C, Huffman D. Absorption and scattering of light by small particles. Weinheim, Germany: Wiley (1983).

CrossRef Full Text | Google Scholar

32. Kuwata H, Tamaru H, Esumi K, Miyano K. Resonant light scattering from metal nanoparticles: practical analysis beyond Rayleigh approximation. Appl Phys Lett (2003) 83:4625–7. doi:10.1063/1.1630351

CrossRef Full Text | Google Scholar

33. Schebarchov D, Auguié B, Le Ru EC. Simple accurate approximations for the optical properties of metallic nanospheres and nanoshells. Phys Chem Chem Phys (2013) 15:4233–42. doi:10.1039/C3CP44124E |

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Moroz A. Depolarization field of spheroidal particles. J Opt Soc Am B (2009) 26:517–27. doi:10.1364/JOSAB.26.000517

CrossRef Full Text | Google Scholar

35. Kelly KL, Coronado E, Zhao LL, Schatz GC. The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment. J Phys Chem B (2003) 107:668–77. doi:10.1021/jp026731y

CrossRef Full Text | Google Scholar

36. González-Castillo JR, Rodriguez E, Jimenez-Villar E, Rodríguez D, Salomon-García I, de Sá GF, et al. Synthesis of ag@silica nanoparticles by assisted laser ablation. Nanoscale Research Letters (2015) 10:399. doi:10.1186/s11671-015-1105-y |

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ermakov VA, Jimenez-Villar E, Silva Filho JMC, Yassitepe E, Mogili NVV, Iikawa F, et al. Size control of silver-core/silica-shell nanoparticles fabricated by laser-ablation-assisted chemical reduction. Langmuir (2017) 33:2257–62. doi:10.1021/acs.langmuir.6b04308 |

PubMed Abstract | CrossRef Full Text | Google Scholar

38. González-Castillo JR, Rodríguez-González E, Jiménez-Villar E, Cesar CL, Andrade-Arvizu JA. Assisted laser ablation: silver/gold nanostructures coated with silica. Appl Nanosci (2017) 7:597–605. doi:10.1007/s13204-017-0599-2

CrossRef Full Text | Google Scholar

39. Pereira MIA, Pereira G, Monteiro CAP, Geraldes CFGC, Cabral Filho PE, Cesar CL, et al. Hydrophilic quantum dots functionalized with gd(iii)-do3a monoamide chelates as bright and effective t1-weighted bimodal nanoprobes. Sci Rep (2019) 9:2341. doi:10.1038/s41598-019-38772-8 |

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bykov DS, Schmidt OA, Euser TG, Russell PS. Flying particle sensors in hollow-core photonic crystal fibre. Nat Photon (2015) 9:461–5. doi:10.1038/nphoton.2015.94

CrossRef Full Text | Google Scholar