Comment: The effect of Mie resonances on trapping in optical tweezers

Bo Sun and David G. Grier
Department of Physics and Center for Soft Matter Research, New York University, New York, NY 10003

Date: February 20, 2009

Abstract:

Recently, Stilgoe, et al., (Optics Express 16, 15039 (2008)) reported calculations of the force on an optically trapped sphere performed using the ``optical tweezers toolbox''. This software suffers from numerical inaccuracies that lead to qualitative errors in the state diagram for stable trapping, particularly for spheres smaller than the wavelength of light.

An optical tweezer is a single-beam optical trap that can localize objects in three dimensions using forces exerted by a strongly focused beam of light (1). Recently, Stilgoe et al. reported calculations of the axial trapping efficiency of optical tweezers for dielectric spheres, which emphasized the important role of Mie resonances on the stability and strength of optical traps, particularly for high-index particles (2). These calculations are based on the Lorenz-Mie theory of light scattering by spheres (3), and so should provide exact results for all sphere sizes and compositions. This approach also is numerically efficient, and can be generalized to compute the forces and torques due to arbitrary light fields (4) on objects of arbitrary shape (6,5). The calculations reported in Ref. (2) were performed with the ``optical tweezers toolbox'' software suite (7) described in Ref. (5), which is formally correct, but relies on numerically inaccurate Hankel functions provided by Matlab (versions R2007a through R2008b, The MathWorks, Natick, MA). Consequently, some of the results reported in Ref. (2) are incorrect, particularly for very small spheres.

Figure 1: State diagram for dielectric spheres in an optical tweezer of numerical aperture 1.3 as a function of the sphere's radius and refractive index. Colors in the main diagram correspond to the maximum axial trapping efficiency at each set of conditions. White regions denote conditions for which particles are not stably trapped. Images (a) and (b) show $\left\vert \vec {E}_i(\vec {r}) \right\vert ^2$ on the surface of the sphere for the particular conditions indicated on the diagram. (a) corresponds to an unstable point. The sphere at (b) is stably trapped.

The incident light field, $\vec{E}_0$ , the scattered wave, $\vec{E}_{\text{s}}$ , and the field within the sphere, $\vec{E}_{\text{i}}$ , are all described in Lorenz-Mie theory as expansions in vector spherical harmonics, (3)

$$\vec{E}_{0}(\vec{r}) = \sum_{n = 1}^\infty \sum_{m=-n}^n a_{nm} \vec{M}^{(2)}_{nm}(k\vec{r}) + b_{nm} \vec{N}^{(2)}_{nm}(k\vec{r}),$$ (1)

$$\vec{E}_{\text{s}}(\vec{r}) = \sum_{n=1}^\infty \sum_{m=-n}^n p_{nm} \vec{M}^{(1)}_{nm}(k\vec{r}) + q_{nm} \vec{N}^{(1)}_{nm}(k\vec{r}),$$ (2)

$$\vec{E}_{\text{i}}(\vec{r}) = \sum_{n=1}^\infty \sum_{m=-n}^n c_{nm} \left( \vec{M}^{(1)}_{nm}(k\vec{r}) + \vec{M}^{(2)}_{nm}(k\vec{r}) \right)$$

$$+ d_{nm} \left( \vec{N}^{(1)}_{nm}(k\vec{r}) + \vec{N}^{(2)}_{nm}(k\vec{r}) \right),$$ (3)

where $k = 2\pi n_m /\lambda$ is the wavenumber of light in a medium of refractive index $n_m$ .

The incident field's coefficients, $p_{nm}$ and $q_{nm}$ , were computed in Ref. (2) through the over-determined point-matching method (8) but can be derived equally well with the Debye-Wolf path integral formalism (4). The scattering coefficients, $a_{nm}$ and $b_{nm}$ , depend on the on the size and composition of the sphere, and on the structure of the light field illuminating it. The interior field's coefficients, $c_{nm}$ and $d_{nm}$ , similarly are computed to satisfy continuity conditions at the sphere's surface (6):

$$c_{nm} = \frac{ i \, m_p \, a_{mn}}{ \xi^\prime_n(x) \, \psi_n(m_px) - m_p \xi_n(x) \, \psi^\prime_n(m_px)}$$ (4)

$$d_{nm} = \frac{ i \, m_p \, b_{mn}}{ m_p \xi^\prime_n(x) \, \psi_n(m_px) - \xi_n(x) \, \psi^\prime_n(m_px)},$$ (5)

where $x = ka$ , $m_p = n_p/n_m$ is the particle's relative refractive index, and where $\psi_n(x)$ and $\xi_n(x)$ are Ricatti-Bessel functions.

Both $\vec{E}_0(\vec{r})$ and $\vec{E}_{\text{s}}(\vec{r})$ contribute to the Maxwell stress tensor on the sphere's surface, whose integral yields the force. References (2) and (5) follow Crichton and Marston (9) by performing the integral analytically. The result for the axial trapping efficiency in an optical trap of total power $P = \sum_{n=1}^\infty \sum_{m=-n}^n\left\vert a_{nm} \right\vert ^2 + \left\vert b_{nm} \right\vert ^2$ is then (9)

$$Q = \frac{2}{P} \sum_{n=1}^\infty \frac{1}{n+1} \sum_{m = -n}^n \... ...q_{nm} q_{n+1,m}^\ast - a_{nm} a_{n+1,m}^\ast - b_{nm} b_{n+1,m}^\ast \right\}.$$ (6)

Equation (10) of Ref. (5) misstates this result, inconsistently taking the real part of the second term. The distributed software, however, is formally correct (7).

The central results of Ref. (2) rely on numerical convergence of the sums in Eq. (6). Unfortunately, the scattering coefficients, $a_n$ and $b_n$ , are notoriously difficult to compute accurately (6,3) because the Ricatti-Bessel functions in which they are expressed have numerically unstable recurrence relations. The expansions' accuracy can be assessed by checking the fields' continuity on the sphere's surface using Eqs. (3) through (5). In our implementation, the relative error in the computed fields is smaller than $10^{-5}$ so that the relative error in the force is less than $10^{-10}$ .

Figure 1 shows the maximum axial trapping efficiency, $Q_{\text{max}}$ as a function of the sphere's radius, $a$ and relative refractive index, $m$ , in an optical tweezer of numerical aperture 1.2. These conditions correspond to those of Fig. 2(b) of Ref. (2). The principal qualitative difference is that Fig. 1 predicts stable trapping for spheres smaller than $a = 0.2\lambda$ , while Ref. (2) does not. The present results are consistent with an analytic computation of the force on a Rayleigh particle in a Gaussian trap (10), and thus would appear to be correct. They also agree quantitatively with results obtained by numerically integrating the Maxwell stress tensor over the surface of the sphere (4), which confirms the correct implementation of Eq. (6). Although these results suggest that even vanishing small spheres can be trapped, the maximum trapping force is limited in practice by saturation of the sphere's polarization, a nonlinear effect not considered in the linear light-scattering formalism discussed here and in Ref. (2).

Explicitly computing the surface fields, such as the examples presented in Fig. 1(a) through (e), reveals smooth symmetric surface field distributions for stably trapped spheres and smooth but asymmetric distributions for unstable conditions. The results in Fig. 1 therefore appear to be free from obvious numerical artifacts.

Major qualitative and quantitative discrepancies we observe with the results of Ref. (2) can be attributed to the computational toolbox's use of the standard Matlab routine besselh(), which implements the Hankel function $H_n^{(2)}(x)$ that is used to compute the spherical Hankel function $h_n^{(2)}(x) = \sqrt{\pi/(2x)} H_{n+\frac{1}{2}}^{(2)}(x) = \xi_n(x)/x$ . Taking $x = 0.118$ as a test case, accurate calculation of $Q_{\text{max}}$ requires the sum in Eq. (6) to be computed to at least $n = 5$ (6,3), with $h_5^{(2)}(0.118) = 2.1996469 \times 10^{-9} - i \, 3.5032872 \times 10^{8}$ . The corresponding result obtained with Matlab is $2.7184 \times 10^{-8} - i \, 3.5033 \times 10^{8}$ , whose real component is in error by a factor of ten. Such numerical errors diminish with increasing $x$ and are negligible for $x > 1$ . Other quantitative differences between Fig. 1 and Fig. 2(b) of Ref. (2) are likely due to slight differences in the description of the incoming wave.

This work was supported in part by the NSF through Grant Number DMR-0606415 and in part by a grant from the Keck Foundation. B.S. acknowledges support from a Kessler Family Foundation Fellowship.