Tutorial on the analytical calculation of optical forces on spherical particles in optical tweezers

TL;DR

This paper provides a comprehensive tutorial on calculating optical forces on spherical particles in optical tweezers using an electromagnetic model that accounts for size, shape, and multipole effects, improving force calibration accuracy.

Contribution

It introduces a complete electromagnetic model for optical forces on spheres, surpassing previous approximations and including effects like resonances and multipole interactions.

Findings

Model accounts for morphology-dependent resonances

Size dependence for large spheres included

Multipole effects for small particles considered

Abstract

Arthur Ashkin was awarded the 2018 Nobel prize in physics for the invention of optical tweezers. Since the first publication in 1986 Optical Tweezers have been used as a tool to measure forces and rheological properties of microscopic systems. For the calibration of these measurements, the knowledge of the forces is fundamental. However, it is still common to deduce the optical forces from assumptions based on the particle size with respect to the trapping laser wavelength. This show the necessity to develop a complete and accurate electromagnetic model that does not depend on early approximations of the force model. Furthermore, the model we developed has several advantages, such as morphology--dependent resonances, size dependence for large spheres, and multipole effects from smaller particles, just to name a few. In this tutorial, we review and discuss the physical modeling of optical forces in optical tweezers, which are the resultant forces exerted by a trapping beam on a sphere of any size and composition.

Figures

Equations

$${\mathbf{E}}_{\text{inc}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(G_{lm}^{TE}{\mathbf{M}}_{lm}^{(1)}(\mathbf{r})+G_{lm}^{TM}{\mathbf{N}}_{lm}^{(1)}(\mathbf{r})\right),$$

$$Z{\mathbf{H}}_{\text{inc}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(G_{lm}^{TM}{\mathbf{M}}_{lm}^{(1)}(\mathbf{r})-G_{lm}^{TE}{\mathbf{N}}_{lm}^{(1)}(\mathbf{r})\right),$$

$${\mathbf{E}}_{\text{sca}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(b_{lm}{\mathbf{M}}_{lm}^{(3)}(\mathbf{r})+a_{lm}{\mathbf{N}}_{lm}^{(3)}(\mathbf{r})\right),$$

$$Z{\mathbf{H}}_{\text{sca}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(a_{lm}{\mathbf{M}}_{lm}^{(3)}(\mathbf{r})-b_{lm}{\mathbf{N}}_{lm}^{(3)}(\mathbf{r})\right),$$

$${\mathbf{E}}_{\text{int}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(d_{lm}{\mathbf{M}}_{lm}^{(1)}(\mathbf{r})+c_{lm}{\mathbf{N}}_{lm}^{(1)}(\mathbf{r})\right),$$

$$Z{\mathbf{H}}_{\text{int}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(c_{lm}{\mathbf{M}}_{lm}^{(1)}(\mathbf{r})-d_{lm}{\mathbf{N}}_{lm}^{(1)}(\mathbf{r})\right),$$

$${\left[{\mathbf{E}}_{\text{inc}}(\mathbf{r})+{\mathbf{E}}_{\text{sca}}(\mathbf{r})-{\mathbf{E}}_{\text{int}}(\mathbf{r})\right]\times \hat{\mathbf{r}}|}_{r=a}=0,$$

$${\left[{\mathbf{H}}_{\text{inc}}(\mathbf{r})+{\mathbf{H}}_{\text{sca}}(\mathbf{r})-{\mathbf{H}}_{\text{int}}(\mathbf{r})\right]\times \hat{\mathbf{r}}|}_{r=a}=0.$$

$$x=k_{2}a=\frac{2\pi n_{2}a}{\lambda },Mx=\frac{n_1}{n_2}{k}_{2}a=k_{1}a,$$

$$\sum _{l=1}^{\infty }\sum _{m=-l}^l\left(G_{lm}^{TE}{j}_l(x)+b_{lm}{h}_l^{(1)}(x)-d_{lm}{j}_l(Mx)\right){\mathbf{X}}_{lm}(\hat{\mathbf{r}})\times \hat{\mathbf{r}}=0.$$

$$G_{lm}^{TE}{j}_l(x)+b_{lm}{h}_l^{(1)}(x)-d_{lm}{j}_l(Mx)=0.$$

$$G_{lm}^{TM}{j}_l(x)+a_{lm}{h}_l^{(1)}(x)-\frac{Z_1}{Z_2}{c}_{lm}{j}_l(Mx)=0.$$

$${\mathbf{N}}_{lm}(\mathbf{r})\times \hat{\mathbf{r}}=\frac{i}{k}\left(\nabla \times (z_l(kr){\mathbf{X}}_{lm}(\hat{\mathbf{r}}))\right)\times \hat{\mathbf{r}}=\frac{i}{kr}\frac{\partial (kr{z}_l(kr))}{\partial kr}.$$

$$G_{lm}^{TM}\frac{\partial (k_{1}r{j}_l(k_{1}r))}{\partial k_{1}r}+a_{lm}\frac{\partial (k_{1}r{h}_l^{(1)})}{\partial k_{1}r}-c_{lm}\frac{1}{M}\frac{\partial (k_{2}r{j}_l(k_{2}r))}{\partial k_{2}r}=0.$$

$$G_{lm}^{TE}\frac{\partial (k_{1}r{j}_l(k_{1}r))}{\partial k_{1}r}+b_{lm}\frac{\partial (k_{1}r{h}_l^{(1)})}{\partial k_{1}r}-\frac{Z_1}{Z_2}{d}_{lm}\frac{1}{M}\frac{\partial (k_{2}r{j}_l(k_{2}r))}{\partial k_{2}r}=0.$$

$${\psi }_l(y)=y{j}_l(y),{\xi }_l(y)=y{h}_l^{(1)}(y).$$

$$G_{lm}^{TM}{\psi }_l^{\prime }(x)+a_{lm}{\xi }_l^{\prime }(x)-c_{lm}\frac{1}{M}{\psi }_l^{\prime }(Mx)=0,$$

$$G_{lm}^{TE}{\psi }_l^{\prime }(x)+b_{lm}{\xi }_l^{\prime }(x)-\frac{Z_1}{Z_2}\frac{1}{M}{d}_{lm}{\psi }_l^{\prime }(Mx)=0.$$

$$\frac{Z_1}{Z_2}\frac{1}{M}=\sqrt{\frac{{\mu }_1{ϵ}_2}{{ϵ}_1{\mu }_2}}\frac{n_1}{n_2}=\sqrt{\frac{{\mu }_1{ϵ}_2}{{ϵ}_1{\mu }_2}}\sqrt{\frac{{\mu }_1{ϵ}_1}{{\mu }_2{ϵ}_2}}\approx 1.$$

$$G_{lm}^{TE}{\psi }_l(x)+b_{lm}{\xi }_l(x)-\frac{1}{M}{d}_{lm}{\psi }_l(Mx)$$

$$G_{lm}^{TM}{\psi }_l(x)+a_{lm}{\xi }_l(x)-c_{lm}{\psi }_l(Mx)$$

$$G_{lm}^{TM}{\psi }_l^{\prime }(x)+a_{lm}{\xi }_l^{\prime }(x)-\frac{1}{M}{c}_{lm}{\psi }_l^{\prime }(Mx)$$

$$G_{lm}^{TE}{\psi }_l^{\prime }(x)+b_{lm}{\xi }_l^{\prime }(x)-d_{lm}{\psi }_l^{\prime }(Mx)$$

$$a_{lm}=-G_{lm}^{TM}\left[\frac{M{\psi }_l^{\prime }(x){\psi }_l(Mx)-{\psi }_l(x){\psi }_l^{\prime }(Mx)}{M{\xi }_l^{\prime }(x){\psi }_l(Mx)-{\xi }_l(x){\psi }_l^{\prime }(Mx)}\right]=-G_{lm}^{TM}{a}_l.$$

$$b_{lm}=-G_{lm}^{TE}\left[\frac{M{\psi }_l(x){\psi }_l^{\prime }(Mx)-{\psi }_l^{\prime }(x){\psi }_l(Mx)}{M{\xi }_l(x){\psi }_l^{\prime }(Mx)-{\xi }_l^{\prime }(x){\psi }_l(Mx)}\right]=-G_{lm}^{TE}{b}_l,$$

$$a_l$$

$$b_l$$

$$\int Y_{lm}^{\ast }{Y}_{pq}d\Omega ={\delta }_{lp}{\delta }_{mq},$$

$$Y_{lm}=\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}{P}_l^m(\cos (\theta ))e^{im\varphi }.$$

$$\mathbf{r}\cdot {\mathbf{E}}_{\text{inc}}(\mathbf{r})=E_0\sum _{l=1}^{\infty }\sum _{m=-l}^l{G}_{lm}^{TM}\mathbf{r}\cdot {\mathbf{N}}_{lm}^{(1)}(\mathbf{r}).$$

$$\mathbf{r}\cdot {\mathbf{N}}_{lm}^{(1)}(\mathbf{r})=-\frac{j_l(kr)}{k}\sqrt{l(l+1)}{Y}_{lm}.$$

$$\mathbf{r}\cdot {\mathbf{E}}_{\text{inc}}(\mathbf{r})=-\frac{E_0}{k}\sum _{l=1}^{\infty }\sum _{m=-l}^l{G}_{lm}^{TM}{j}_l(kr)\sqrt{l(l+1)}{Y}_{lm}.$$

Full text

††journal: josab††articletype: Research Article

1 Optical trapping

Optical trapping is a modern tool that employs a structured laser beam to manipulate and investigate cellular-scale physio-chemical phenomena. A few of the advantages of the trapping technique are that it exerts no physical contact and exhibits a very low stiffness, which makes it an ideal tool for soft-matter research. Moreover, the trapping beam can be split or multiplexed into several traps. During the last couple of decades, optical trapping has been greatly perfected owing to instrumental development and calibration techniques. Among the different optical trapping techniques, the single beam gradient optical trap is the prevailing one. This trap is also referred to as optical tweezers, a term coined by the technique´s pioneer Arthur Ashkin [1, 2], and acknowledged with a shared Noble prize in physics in 2018. Ashkin’s book [3], a collection of early articles with commentaries, illustrates the early years of optical tweezers research, along with some of the related modeling briefly mentioned in this tutorial.

The physical principle behind optical trapping lies in, the momentum of the electromagnetic field of the trapping beam, as stated by Maxwell. When this beam is incident upon an object, it is scattered elastically upon their encounter, resulting in a net change to the rate of linear momentum. Because of the conservation of linear momentum, a force is exerted on the object by the light beam. To obtain an intuitive idea of the force that can be exerted by light, the light’s momenta is related to de Broglie relation from modern physics, which indicates a minuscule momentum per photon, resulting in a minuscule force, comparable to the weight of cellular-sized objects. For example, a 5 $\mu$m water droplet weight just as much as the force imparted on a mirror by a 1mW laser pointer.

To achieve a significant change in the linear momentum, a highly focused beam is required. In standard optical tweezers, this trapping beam is a single highly focused Gaussian profile laser beam. The trapped object is generally a dielectric microsphere, for which forces can be calibrated and employed as a force transducer to determine optical forces in other non-spherical samples within the same medium.

1.1 Optical force models

Early theoretical modeling for optical forces conducted as early as 1909 by Debye, showing that the force on a sphere can be described in terms of extinction and scattering cross sections [4, 5]. Debye considered a linearly polarized plane wave incident on a homogeneous spherical particle.

To theoretically determine these optical forces, initial assumptions are set based on the microsphere size with respect to the beam wavelength. For a small object that is much smaller than the beam wavelength, the object behaves like an electric dipole, in agreement with Rayleigh scattering. Hence, a simple description is established for optical forces acting on an electromagnetic dipole, and this is known as the Rayleigh regime [6]. However, for objects greater than the wavelength, the trapping beam is treated as rays. The optical force is thus described in terms of geometrical optics, considering multiple reflections of the rays in the trapped microsphere to determine the net change in momentum [7, 8]. However, in the case of most experiments, the trapped object diameter is within the order of the laser wavelength, and the above mentioned approximations fail. A complete classical electromagnetic (EM) wave theory is needed to model optical forces in this intermediary regime, known as the Mie regime, although the Mie theory is applicable to spheres of any size. This analytical EM approach has been avoided mainly owing to its troublesome calculations; however, after some algebraic derivations, simple expressions for accurate optical forces, valid for any object size, can be obtained.

The complete EM approach allows to determine the morphology–dependent resonance, or whispering gallery modes, which have been observed in optical traps and are absent in the approximate size regime descriptions [9, 10, 11]. The early studies that considered the complete electromagnetic approach are those of Mie [12] and Lorenz [13]. Historical reviews of the initial developments and contributions are provided by Logan [14] and Kerker [5]. The Lorenz-Mie theory is an analytical method designed for plane-wave scattering by a spherical particle, which is not applicable to optical tweezers that require highly focused beams. Therefore, a generalization for the scattering of generic beams by Mie sized particles led to the development of the generalized Lorenz-Mie theory (GLMT)[15, 16, 17]. In the GLMT framework, the incident and scattered beams are expressed in terms of vector spherical wavefunctions, whose expansion coefficients describe the applied beam configuration.

Once the vector scattered fields are known, we can, therefore, proceed with the Maxwell stress tensor (MST) approach to our problem and determine the net force or torque the beam exerts on the microsphere. This approach was first performed by Kim and Lee [18], where they presented a vector theory for optical forces acting on a sphere by a Gaussian beam described from a complex-source-point method. The treatment relied on expanding the fields in vector spherical harmonics, much like GLMT. Subsequently and in a similar manner, Marston and Crichton [19], presented an ab initio calculation of the radiation torque on an isotropic sphere illuminated by a circularly-polarized plane electromagnetic wave. The same approach was soon thereafter extended to optical forces and torques by Barton et. al [20], considering a focused beam model from Davis [21, 22]. The Davis beam model depicts a paraxial Gaussian beam with higher order corrections.

From the initial development of GLMT, most applications required knowledge of the scattered field or beam shape coefficients, which are none other than the multipole expansion coefficients. Several strategies have been developed for the determination of the beam shape coefficients (BSCs) in the most efficient and/or accurate way.

Apart from the several ways to determine the BSCs, at the heart of this development, the quest of the appropriate incident beam model was also highlighted. The suggested and employed beam models were not a solution to Maxwell’s equation, and most of them were approximations that broke down for a highly focused beam.

Many researchers contributed to the understanding of optical forces by ways of theory and computational modeling [23, 24, 25, 26, 27, 28, 29, 17, 30]. In this tutorial, we wish to present what we believe is a didactic approach towards an analytical solution, considering a physical and realistic trapping beam along with the exact determination of its beam shape coefficients. With the MST approach, an exact model for optical forces in a typical optical tweezers is obtained [11]. Therefore, we merge the main contributions from some of the studies highlighted in this work, into what we believe is a robust and rigorous treatment for optical forces in optical tweezers. Similar approaches using BSC calculation via plane wave spectrum or for Gaussian beams independent of the radial coordinate, have been published [31, 32, 33, 34], but were unaware of the proposed method in Ref. [11],

A standard and modern notation for the equations was adopted throughout this work. The tutorial is divided into four parts. The first part introduces the EM fields along with the multipole coefficients. The second part explores how the BSCs are determined, and this is illustrated for the simple case of a plane wave. The third part presents the highly focused beam description along with the BSCs. Finally, the optical force through the Maxwell stress tensor is presented in terms of the BSCs. We hope to have addressed all the important points and developments described above.

2 Modeling

Hence, how do we start modeling the optical forces? Primarily, we rely on the knowledge of the exact electromagnetic (EM) fields, which travel as a wave with linear momentum. Because the trapped object is a microsphere, it seems reasonable to express the EM fields in terms of spherical harmonic functions, similarly to the standard textbook approach [35]. The encountered problem is essentially the elastic scattering of monochromatic coherent light by the homogeneous and isotropic sphere of radius ($a$).

2.1 Description of the EM fields

Let us describe all the required electromagnetic fields, incident (inc), scattered (sca), and internal (int) as time-harmonics in terms of spherical vector expansion, also known as multipole expansion, as our scattering particle has spherical symmetry [36]. The time dependence $e^{-i\omega t}$ is omitted for clarity. Our coordinate system has its origin at the center of our spherical scatterer, which scatters the incoming light. A similar approach with other notations is present in several textbooks [5, 37].

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $l$ and $m$ are the radial and azimuthal mode index, respectively. The electric field dimension constant $E_{0}$ is used to make the expansion coefficients non–dimensional. The mode amplitudes, $a_{lm},b_{lm},c_{lm},d_{lm},G_{lm}^{TM},$ and $G_{lm}^{TE}$, are also called expansion coefficients or multipole coefficients, however, the last pair is special, and we will refer to them as beam shape coefficients (BSCs). The vector spherical harmonics are, $k\mathbf{N}_{lm}(\mathbf{r})=i\nabla\times\mathbf{M}_{lm}(\mathbf{r})$, $\mathbf{M}_{lm}(\mathbf{r})=z_{l}(kr)\mathbf{X}_{lm}(\mathbf{\hat{r}})$, and $\mathbf{X}(\mathbf{\hat{r}})=\mathbf{L}Y_{lm}(\mathbf{\hat{r}})/\sqrt{l(l+1)}$ [38], where $z_{l}(kr)$ is the appropriate spherical Bessel function, $Y_{lm}(\mathbf{\hat{r}})$ are the scalar spherical harmonics, $Z=\sqrt{\mu/\varepsilon}$, $k=\omega\sqrt{\mu\varepsilon}$, and $\mathbf{L}=-i\mathbf{r}\times\nabla$ is the angular momentum operator in direct space. The above vector functions are solutions to the vector wave equation, first proposed by Hansen [39], and any field can be expanded in these series given the appropriate coefficients.

The spatial function $z^{(j)}_{l}(kr)$ is a basic solution to the spherical Bessel equation, which for $j=1,2,3,4$ corresponds to the ordinary spherical Bessel, spherical Neumann, and spherical Hankel functions of the first and second kind, respectively. The choice of the function depends on the particular spatial bounds, either for the fields inside or outside the scatterer. Because the fields must be regular at the origin, only $z_{l}^{(1)}(kr)=j_{l}(k_{1}r)$ can appear at the interior radial function, and also for the incident fields that exist in the origin, independent of the scatterer. While in the region outside the scatterer, we can write in terms of spherical Hankel functions of the outgoing type, based on the physical ground of the scattered wave, $z_{l}^{(3)}(kr)=h_{l}^{(1)}(k_{1}r)$ .

2.2 Scattering coefficients

To solve the electromagnetic boundary problem on the scatterer surface, the continuity of the tangential components of the electric and magnetic fields is required. The following equations need to be satisfied,

[TABLE]

[TABLE]

To solve each of these tangential components, let us introduce two new abbreviations,

[TABLE]

where $x$ is the size parameter, $M$ is the relative refractive index, $n_{1}$ and $n_{2}$ are the refractive index of the particle and the surrounding medium, respectively. Note that the refractive index of the particle $n_{1}$ is allowed to be complex, to treat the effects of absorption such as the ones in metallic particles. For the electric TE modes, the solution is

[TABLE]

This is valid for all angles, therefore

[TABLE]

Analogously, for the magnetic TE modes, the condition valid for all angles is as follows

[TABLE]

The terms proportional to $\mathbf{N}_{lm}(\mathbf{r})$ simplify to the following

[TABLE]

For the TM modes, we obtain

[TABLE]

Moreover, for the magnetic equation, we have

[TABLE]

To further simplify the expressions above, we introduce the Riccati–Bessel functions,

[TABLE]

Which leads to

[TABLE]

and

[TABLE]

Considering the relations between the EM impedance and the waves at optical frequencies $\mu_{1}\approx\mu_{2}$,

[TABLE]

The resulting four equations relating all multipole coefficients are:

[TABLE]

Because we are interested in the scattered fields, we can eliminate $c_{lm}$ and $d_{lm}$. First, we isolate $c_{lm}$ in Eq.(21) and substitute in Eq.(22).

[TABLE]

Secondly, we isolate $d_{lm}$ in Eq.(23) and substitute in Eq.(20).

[TABLE]

where the Mie scattering coefficients are thus defined as in Ref.[37],

[TABLE]

A complete and efficient description of the scattered EM wave can be obtained if information about the scatterer and the incident field are available. The GLMT framework addresses the fields in all space (near and far fields), allowing for a complicated beam that satisfies Maxwell’s equations. In our case, the scatterer is a homogeneous and isotropic dielectric microsphere, and it is completely characterized by the Mie scattering coefficients ($a_{l}$ and $b_{l}$), through the knowledge of the microsphere radius ($a$), the refractive index of the microsphere ($n_{2}$), the refractive index of the surrounding medium ($n_{1}$), and the wavelength of the incident wave ($\lambda$). For the incident beam, this information is present in the beam shaped coefficients ($G_{lm}^{TM}$ and $G_{lm}^{TE}$), which have been studied extensively by Gousbet, Grehan, and Lock [40, 17] (and references therein). The method by which we determine the beam shaped coefficients is presented in the following section.

3 Beam shape coefficients

The description of BSCs appears in the GLMT framework, which was initially pioneered when studying plasma velocity by laser Doppler velocimetry [17]. The results for the scattering of any beam are known once the beam shape coefficient is determined. In fact, GLMT can be applied to analytically solve scattering from other objects besides the homogeneous sphere, including the stratified sphere, spheroids, and infinite cylinders of any size. Moreover, any type of beam can be used, the solution not being restricted to Gaussian beams. The main problem in many studies and the subject of much debate is the calculation of the BSCs, mainly due to the non–cancellation of the radial spherical Bessel functions in both sides of the expansion. Primarily, they are defined as follows.

In Equations (1) to (6) we described all electromagnetic fields in terms of spherical vector expansion. In the previous section, we explicitly determined the complex scattering coefficients in terms of the beam shape coefficients (BSCs) $G_{lm}^{TM}$ and $G_{lm}^{TE}$. We can, therefore, express any incident EM field by determining these BSCs. Note that $\mathbf{M}_{lm}(\mathbf{r})$ is perpendicular to the radial vector, $\mathbf{\hat{r}}$, we can use this and the normalization relation of spherical harmonics;

[TABLE]

where

[TABLE]

To determine the BSCs $G_{lm}^{TM}$ and $G_{lm}^{TE}$, we first multiply by the radial vector $\mathbf{r}$;

[TABLE]

The last vector product in the equation above is determined by expanding in terms of spherical vector harmonics, and later in terms of spherical harmonics,

[TABLE]

Therefore, Eq.(30) becomes

[TABLE]

To eliminate the double sum, we will multiply both sides by $Y_{pq}^{*}$ and integrate over the solid angle as follows

[TABLE]

Using the normalization property of Eq.(28), we obtain

[TABLE]

or

[TABLE]

Analogous for the magnetic field, the corresponding BSC is given by

[TABLE]

The Eqs.(35) and (36) above define the beam shape coefficients for the incident trapping beam from the GLMT framework. The right–hand side of these equations should contain a spherical Bessel function from the integration to cancel with the one on the left–hand side. The major shortcoming of the other methods is that this radial dependent term has never been explicitly determined to produce BSC independent of the radial coordinates.

The radial component of the electric and magnetic fields of the incident beam needs to be determined. Once the exact EM is known, or described an analytical expression, we are able to recover an exact BSCs. However, if the incident EM fields is only approximated, we no longer retain an exact BSC. This is the case for the Davis beam, which results in a BSCs that depends on the radial coordinate ($r$). By definition, the expansion coefficient should not depend on the radial coordinate, hence Barton suggested to determine the BSCs at the particle surface ($r=a$) [20]. However, what would happen if $ka$ is a root of the spherical Bessel function ($j_{l}(ka)$)? Note that there is an infinite sum in the order $l$ and infinite roots for each order.

3.1 BSC for a plane wave

The simplest incident field is a plane wave. Here, we determine the BSCs of the plane wave by direct integration. Consider an incident coherent monochromatic plane wave propagating in non-absorbing homogeneous media in the z–direction, which is linearly polarized in the x–axis,

[TABLE]

In spherical coordinates, the radial component is

[TABLE]

We can now calculate the plane wave BSCs from the definition of Eq.(35),

[TABLE]

The solution for the integral in $\phi$ is,

[TABLE]

Therefore, for plane waves only BSCs with $m=\pm 1$ remains. Let us determine the BSC for $m=+1$.

[TABLE]

The $\theta$–integral can be solved using the integral representation of the spherical Bessel function,

[TABLE]

Eq.(41) thus becomes

[TABLE]

For $m=-1$, we have

[TABLE]

Similarly, for the other BSC,

[TABLE]

The $\phi$ integral now becomes

[TABLE]

For $m=\pm 1$,

[TABLE]

We can verify that these are the true BSC, by inserting them into the radial component of the incident electric field Eq.(1).

[TABLE]

The last summation in Eq.(48) is a well–known relation, known as the plane wave expansion or Rayleigh equation, which describes the plane wave as a sum of spherical waves. As a result, we recover the expected electric field of Eq.(38). Having explained how to obtain the plane wave BSC, let us turn to the beam of interest for optical tweezers, namely the highly focused Gaussian beam.

4 Highly focused Gaussian beam

Another great challenge for a proper beam description is using an expression for the EM field that satisfies exactly Maxwell’s Equation. The beam representation must represent a solution of the vector Helmholtz equation, and the paraxial Gaussian beam is not a solution. The use of these pseudo-fields [41], or non-Maxwellian beams, would result in an approximate beam description leading to a pair of BSCs (Eqs.(35) and (36)) that depend on the radial coordinate [42].

Early studies modeled a focused Gaussian beam as a paraxial beam perturbed up to the fifth-order through a series expansion in terms of a $s-$parameter ($s=1/(kw_{0})$, where $w_{0}$ is the beam waist at the focus) [22]. This Davis-Barton fifth-order approximation still does not solve Maxwell’s equation exactly and fails to converge to reproduced experimental results under highly focused conditions ($s\rightarrow 1$) as is the case in optical trapping [43, 44]. The main limitation of beam models is that they do not satisfy exactly Maxwell’s equation, and do not represent a good approximation of real beam when the spot size is too small. Most paraxial beams used in focusing are also non–physical because they do not exhibit a change in spot size due to polarization alone [45]. Moreover, the diffraction effects due to the truncated Gaussian beam by the finite objective lens and the polarization mixing by the high numerical aperture (NA) objective are not correctly represented near the focus of a high-NA objective [24, 46].

The suggested most efficient way to overcome this issue of an accurate beam model for optical force prediction was to treat the high-NA focusing as a vector diffraction problem [43].

At high-NA focusing, the paraxial scalar theory does not take into account the vector nature of the beam, it does not justify small-angle approximations, and it does not consider the apodization effects of optical systems [47, 48]. This has been thoroughly tested, proven and compared in numerous experiments, such as the 4Pi confocal fluorescence microscope [49], which led to the development of super–resolved fluorescence microscopy. A tutorial for diffraction calculation in optical systems has recently been published [50], discussing diffraction integrals and high NA focusing.

The solution for the highly focused beam that represents the vector property of the incident beam is the vectorial Debye–Wolf diffraction integral method, proposed by Richards and Wolf [51]. This model describes the paraxial Gaussian beam being focused by a highly convergent aplanatic lens as a sum of plane waves arriving from all the angles spanned by the lens´s numerical aperture as an angular spectrum of plane waves [52]. Recently, it was highlighted that this approach to vector diffraction dates to a much earlier publication, in 1920 by Ignatovsky [53, 54].

4.1 Angular Spectrum representation of focused Gaussian beams

The angular spectrum representation depicts a mathematical instrument to calculate the electromagnetic fields near a highly focused beam (by a lens of focal length $f$). It comprises the decomposition of the EM field into an infinite series expansion in terms of Fourier components. With respect to optical tweezers, we are concerned with the propagating part of the angular spectrum, which can be solved with the method of stationary phase [55].

The resulting electric field in the vicinity of the focus ($\mathbf{E_{inc}}(\rho,\varphi,z)$), as lens focus is chosen as the origin of the coordinate system, with respect to the field on a reference sphere ($\mathbf{E_{\infty}}(\theta,\phi)$) is given in cylindrical coordinates as follows:

[TABLE]

where $\theta_{\text{NA}}=\sin^{-1}(\text{NA}/n_{a})$ is the maximum acceptance angle given by the NA of the objective lens. The field on the reference sphere ($\mathbf{E_{\infty}}$) is itself determined from the paraxial Gaussian $TEM_{00}$ beam as

[TABLE]

where $n_{b}$ and $n_{a}$ refer to the refractive index of the medium before and after the focusing lens, respectively. The electric field in the focal region, Eq.(49), is described in cylindrical coordinates and completely satisfies Maxwell’s equation. The angular spectrum of a Gaussian beam requires a corresponding Gaussian spectrum that includes evanescent components. These non-homogeneous components decay exponentially with distance from the reference sphere. In optical tweezers, the focal region would be far-field from the reference sphere. Therefore, this method allows for a more rigorous description of the focused field in comparison to the previous ones [52].

To determine the BSCs, we need to convert to spherical coordinates. We would also like to displace this beam from its origin (the focal point) to an arbitrary point in space ($\rho_{0}$,$\varphi_{0}$,$z_{0}$). This is achieved with the following substitutions:

[TABLE]

The purpose of shifting the EM field with respect to the sphere origin provides a later advantage in describing the scattering, or the corresponding optical forces, when the focused beam is displaced with respect to the sphere center. In general particle scattering theory, such as the transition matrix (T-matrix) is formulated in the coordinate system located at the particle center [56, 57]. The beam (and it’s BSCs) is described with respect to another origin with the use of the translational–addition theorem [58, 59, 60], which adds another double sum operation that needs to be truncated appropriately. In our case, this is avoided by first generalizing the focused beam description (and its corresponding BSCs) with respect to any arbitrary location.

We now modify the spherical angular coordinates of the focusing system output pupil to avoid confusion with the BSCs integration variables, from ($\theta$,$\phi$) to ($\alpha$,$\beta$). The focused field becomes

[TABLE]

The amplitude of the incoming field, before focusing, has a Gaussian distribution and general polarization:

[TABLE]

where $\omega$ is the beam waist at the objective aperture, and ($p_{x}$,$p_{y}$) are the complex polarization parameters describing the state of polarization of the beam before entering an aplanatic objective lens. For example, linear polarization in the x-axis corresponds to (1,0). Consequently, the magnetic field is given as follows:

[TABLE]

That is, to switch to the BSC which depends on the magnetic field we have to change $p_{x}\rightarrow-p_{y}$ and $p_{y}\rightarrow p_{x}$. Hence, we can determine Eq.(50) with the above incident field, changing to the new $\alpha$ and $\beta$ angle coordinates. For the BSC, we only require the radial component of the incident field, which implies the radial component of the field at the reference sphere,

[TABLE]

Thus, the field at the reference sphere and the radial component at the focus are,

[TABLE]

[TABLE]

These double integrals are addressed in the following section.

4.2 Beam shape coefficients of focused Gaussian beam

We can now determine the BSCs using the definition in Eq.(35), with the radial component of the incident field described in Eq.(59).

[TABLE]

We proceed by solving three of the four integrals above. The first integral to solve is that of $\phi$, where we use the following relation,

[TABLE]

where we abbreviated $y=kr\sin\alpha\sin\theta$. This solves the first of four integrals, and we will further introduce some dummy variables to shorten the algebra.

[TABLE]

which after rearranging term proportional to $\cos\alpha$ becomes

[TABLE]

Both $C$ terms can be combined within the $\theta$ integral, considering it as a result of the derivative with respect to $\alpha$,

[TABLE]

The result for the $\theta$–integral is given in Ref. [61] as follows,

[TABLE]

This reduces the BSC to only two integrals,

[TABLE]

Note that this last $\theta$-integral has actually cancelled the radial dependence from the spherical Bessel function defined in the BSC expression of Eq.(35). The analytical result for this integral is the key to avoid most numerical approximation in literature. This integral has been pointed out as a special case of a Gegenbauer integral [62], even though it is not generally well known and not easily available in integral tables [63]. This result has been useful in several applications other than light scattering, including acoustics [64, 65], ionization rates for orbitals [66], the Casimir interaction [67], atom tunneling [68], field-theory [69], and cosmology [70].

This analytical calculation was not derived for many years, and other strategies have been employed to solve the problem of determining BSCs, such as the use of van de Hulst’s localization principle, building a localized model of the beam. For an incident Gaussian beam on-axis to a sphere, described by the Davis model, the transverse and angular (polar) coordinates are mapped to specific values within the localized approximation as, $kr\rightarrow n+1/2$ and $\theta\rightarrow\pi/2$, respectively, thereby no longer depending on the radial variable, and avoiding the integral over the polar component [71, 72, 73]. A historical overview of the different methods for determining the BSC can be found in Gouesbet’s book [17].

To proceed with the determination of BSCs, we integrate (with respect to $\beta$), as we need to eliminate our dummy indexes that are functions of $\beta$. We introduce two angular functions related to the associated Legendre function, which are common in light scattering [5],

[TABLE]

After some algebra we obtain,

[TABLE]

The integral (with respect to $\beta$) is equal to,

[TABLE]

where $x=k\rho_{0}\sin\alpha$. Solving this last analytical integral, we obtain the BSC for the highly focused incident beam.

[TABLE]

The analogous TE BSC,

[TABLE]

This numerical integration is quite fast once the integrating functions are known. For the calculation of optical forces with the beam positioned in space ($\rho_{0}$,$\phi_{0}$,$z_{0}$) with respect to the sphere center, only the BSCs need to be recalculated, which is an advantage with respect to other methods that do not share this feature. These BSCs can be expressed in a more compact form, uncoupling the azimuthal coordinate ($\phi_{0}$) and polarization ($p_{x}$,$p_{y}$) [74],

[TABLE]

The BSCs $G_{lm}^{TM,TE}$ corresponds to the transverse magnetic (TM) and transverse electric (TE) multipoles, which correspond to the left and –right polarizations of the last vector in Eq.(72), respectively, whose last lines involve a (1$\times$2)(2$\times$2)(2$\times$1) matrix operation. The efficiency of this calculation can be higher by exploiting the symmetry of some special cases, as highlighted in the following sections.

4.3 Case study: coincident microsphere and beam axes

This case is a common condition under trapping equilibrium, resulting in the microsphere axis coincident with the beam axis. To determine these special BSCs, we make, $x\rightarrow 0$, which leads to the following conditions;

[TABLE]

For $m=\pm 1$ we have

[TABLE]

while for $G^{TE}$ we obtain

[TABLE]

that can be summed up as follows:

[TABLE]

where

[TABLE]

5 Optical Forces

The optical force is derived from the conservation law of electromagnetic wave momentum. For a trapped spherical scatterer, the total electromagnetic force may be determined by applying Maxwell’s stress tensor (MST) in the space surrounding the scatterer. A step-by-step application of the MST for an arbitrary incident wave and a dielectric sphere is presented. With the EM field described in terms of spherical harmonics, the observable effects of the force components that an EM wave exerts on the microsphere is given by the time average of the Maxwell stress tensor ($T_{ij}$), which for steady–state optical conditions is given in the traditional Minkowski form, defined below as in Ref.[35],

[TABLE]

where the EM fields are the ones outside the scatterer (incoming and scattered fields). We integrate over a spherical surface surrounding the scatterer, which can be chosen arbitrarily; hence, we choose one whose radius $kr\rightarrow\infty$. Where $n_{j}=\mathbf{\hat{r}}$, and because $dA=r^{2}d\Omega$, only certain terms of the radial function expansion is needed. At infinity, the electromagnetic fields only contain transverse components, and therefore the first two terms of the right-hand side of Eq.(78) do not contribute, resulting in

[TABLE]

The EM fields involve spherical vector terms such as

[TABLE]

Which for $kr\rightarrow\infty$ would result in

[TABLE]

The asymptotic expansion of the spherical radial functions and their derivatives in this limit are

[TABLE]

The total EM field outside the sphere and the sum of the incident and scattered fields are known in terms of the spherical vector harmonic. Inserting these fields into Eq.(79), and taking into account that the incident field contribution is removed as it is independent of the sphere position, the contribution of both fields to the force amounts to zero, i.e., all terms that do not contain $a_{l}$ and $b_{l}$. The vectorial spherical harmonic products can be re-written as follows,

[TABLE]

and

[TABLE]

These simplify the electric and magnetic terms to two components each. In this way, we can sum the contributions from the electric and magnetic terms, resulting in a force

[TABLE]

The force is generally split into components, such that the explicit integrals over the solid angle are given in terms of

[TABLE]

[TABLE]

[TABLE]

[TABLE]

After substituting these integrals in Eq.(88), we obtain after some simplification the transverse and longitudinal force components as follows:

[TABLE]

[TABLE]

where the combination of Mie coefficients has been defined as the following

[TABLE]

The force calculations, outlined step-by-step above, can easily be generalized for optical torques using the same stress tensor and guidelines presented in this tutorial. This methodology was used to calculate the force in an optical force spectroscopy experiment [75]. In Ref. [42], the figures show a comparison between the expansion method and the real beam. This model also explained the results from a double tweezers setup to perform an ultra–sensitive force measurement of the full optical force profile as a function of radial position and wavelength for both transverse electric (TE) and transverse magnetic (TM), modes excitation [11]. Aberrations can also be considered, where the model successfully explained the ripple structure in the optical forces profile due to spherical aberrations due to refractive index mismatching [74].

6 Final Remarks

The infinite summation for the optical force may be of concern. Computationally, we truncate it to a particular value ($N_{max}$), because all other terms have a negligible contribution. This maximum summation term scales with the scatterer size (size parameter, $x$), just as in the determination of the scattering cross section, which uses the Wiscombe criterion [76]. By using a different floating-point variable for the numerical computation, a suitable truncation number can be chosen [77], by taking the closest integer of

[TABLE]

where $10^{-\epsilon}$ is the relative truncation error, and $x$ is the size parameter defined in Eq.(9). This is correct only for the total cross section, as has been emphasized recently [78], or for EM fields within a region of space delimited by $k_{2}r<N_{max}$. However, for specific fields instead of cross sections, a similar rule is applicable in which the maximum number of terms should be greater than the largest $k_{2}r$ in the plot region.

For most applications, when the bead is near its equilibrium position in the optical trap, the trap stiffness ($\mathbf{\kappa}$) may suffice, which can be directly calculated from the optical force from Hooke’s law,

[TABLE]

where $q=x,y,\text{ or }z$, and is evaluated at the equilibrium position ($r_{eq}$). Many refinements have been made to render the approach described here more complete. One can take into account the spherical aberrations produced by the refractive index mismatch at the coverslip interface [79], whose consequences on optical forces are known [80, 81, 82]. The appropriate aberrated beam coefficients can be determined in a similar way [74], making the optical force more realistic when trapping close to an interface.

One of the main efficiencies in determining optical forces with this procedure, are the analytically determined BSCs at an arbitrary position with respect to the microsphere, without the need of introducing computational error from translation operations. The present method allows to obtain the resultant force on a spherical dielectric by displacing the particle with respect to beam focus, and hence a full force profile can be obtained.

The precise determination of the BSC for a highly focused beam has been validated experimentally, making it a robust strategy for BSC determination. One of these validations was the measurement of a complete force profile for microspheres of different sizes, and the effects of coupling the beam polarization to different TM and TE modes separately [11]. Recently, this force profile or force profile mapping measurement and comparison have been performed with an improved setup [83, 84, 85]. The analytical BSCs was also validated independently in a single particle photo-thermal experiment [86].

Another method to analytically determine the BSCs was derived by Moreira et. al using the Fourier transform of the angular momentum operator [87, 88]. In this work, the BSCs are for the first time shown to be independent of the radial coordinate if the beam model exactly satisfies Maxwell’s wave equations. Several, analytical BSC for different beam modes, such as Bessel beams, circular, and rectangular wave-guided beams were likewise derived. This new method also applies to plane waves and highly focused beams.

As a final contribution of this study, the computer codes generated as a result of all of the work presented in this paper are released and are freely available, also allowing for contributions from the community[89].

7 Summary

We presented an overview of the approach to model optical forces occurring in optical tweezers. The optical forces are determined by the Maxwell stress tensor approach and expressed in terms of the beam shape coefficients from the Generalized Lorenz-Mie Theory. A robust beam model is employed for the tightly focused Gaussian beam with a high NA objective lens in terms of the diffraction integrals presented by Richards and Wolf and expressed by the angular spectrum representation. From this beam, the beam shape coefficients are determined by analytically integrating under the solid angle.

This work allows for accurate and precise computational modeling to quantitatively compare the theoretical results with experimental outcomes. Using the analytical expressions presented in this work requires low computational complexity. We hope that this work can be of help to help young researchers in the field and provide them with a more solid background in optical forces theory, stimulating further development for trapping applications.

Funding

This work was partially funded by the São Paulo Research Foundation (FAPESP) (2017/13543-3).