Application of the discrete dipole approximation to optical trapping calculations of inhomogeneous and anisotropic particles

Stephen H. Simpson; Simon Hanna

doi:10.1364/OE.19.016526

1. Introduction

Many techniques are available for computing the forces and torques that arise in optical trapping experiments. Most approaches require a solution to the optical scattering problem, followed by an evaluation of the rate of change of linear or angular optical momentum. These quantities are, in turn, associated with the forces and torques on the particle.

The choice of computational technique depends on the system under study. For particles that are large compared with the optical wavelength, ray tracing may be appropriate [1], while for particles much smaller than the wavelength, the Rayleigh approximation is of use [2]. At intermediate length scales, rigorous solutions of Maxwell’s equations are required. The techniques that have been applied in this regime include the finite element method (FEM) [3], the finite-difference time-domain (FDTD) method [4, 5], the discrete dipole approximation (DDA) [6–11], and T-matrix theory [12–14].

The FEM, FDTD and DDA approaches have the advantage that they can model arbitrary geometries and material parameters without significant modification. However, all three techniques are normally regarded as extremely demanding of computing resources. On the other hand, the T-matrix theory, while extremely powerful, is more restricted in its application owing to its reliance on central expansions of the electric field [15]. There are difficulties when using T-matrix theory to model objects of large aspect ratio, because of the nature of the expansion [15], and arbitrary inhomogeneity and irregularity require some effort to include [16].

In this article we examine the use of an in-house implementation of the full DDA, specifically intended for optical trapping problems, which takes advantage of FFT methods to accelerate the calculation [9]. As its name suggests, the DDA is an approximation – the forces and torques calculated are based on assumptions that have not been rigorously justified (see below, Section 2.2). Therefore, the accuracy and convergence properties of the DDA are assessed by comparing results with exact T-matrix calculations for a range of systems with different symmetries. A similar comparison was presented recently [10], but was restricted to calculations of forces on dielectric spheres. Here we extend the verification of the DDA technique to calculations of forces and torques on lower symmetry particles, including the effects of optical and geometric anisotropy as well as structural inhomogeneity. Thus we consider dielectric spheres and spheroids, birefringent spheres, and dielectric core-shell particles. The results presented include comparative calculations of torques on nonspherical, birefringent and inhomogenous particles, that are not found elsewhere in the literature.

Finally, the generality of the DDA method is demonstrated by performing additional calculations on more general systems of topical interest. We examine the behaviour of the letter-shaped colloidal particles investigated experimentally by Wilking et al. [17], birefringent quartz cylinders that have been used to probe the torsional dynamics of DNA [18, 19], and biphasic Janus particles, which have a variety of potential applications including use as biosensors [20].

2. Methods

Since the DDA and T-matrix methods are being compared, we give a brief description of our implementation of each, followed by a discussion of the generalised stiffness tensor and a short explanation of force and torque densities.

2.1. T-matrix theory

An extensive account of T-matrix theory appears elsewhere [13]. The method involves expansions of incident, scattered and internal fields in an orthogonal basis set, usually the vector spherical wave functions (VSWF). The incident field may be written as:

E^{i n c} = \sum_{n = 1}^{\infty} \sum_{m = - n}^{n} (a_{n, m} Rg M_{n, m} + b_{n, m} Rg N_{n, m})

where RgM _n,m and RgN _n,m are regularised VSWFs. Application of appropriate boundary conditions leads to a linear relationship between the expansion coefficients for the incident and scattered fields, i.e. the T-matrix, thereby solving the scattering problem. For homogeneous and isotropic spheres, the T-matrix is diagonal, with entries given by the Mie coefficients. For non-spherical objects we calculate the T-matrix using the exact, extended boundary condition method (EBCM). Optical anisotropy is incorporated by using the VSWF basis set to construct a new basis in the anisotropic material such that the divergence of the electric displacement vector D vanishes, while simultaneously satisfying an appropriate wave equation [21, 22].

Following solution of the scattering problem, forces are obtained by integration of the Maxwell stress tensor over a closed surface surrounding the scatterer [23]. This integration is performed analytically in the far field, leading to a quadratic sum over the VSWF expansion coefficients. The incident field is provided by a complex source point (CSP) beam [24]. The appropriate expansion coefficients for CSP beams are found by expressing the incident field [9] in spherical polar coordinates and applying inner product formulae [13]:

a_{n, 1} = a_{n, - 1} = b_{n, 1} = - b_{n, - 1} = I_{n, m}

m is restricted to ±1 and I_n,m is given by the following, which is evaluated numerically:

I_{n, m} = i^{n + 1} γ_{n, m} \int_{0}^{π} d θ e^{κ z_{0} cos θ} (1 + cos θ) [P_{n}^{m} (cos θ) + sin θ \frac{d}{d θ} P_{n}^{m} (cos θ)]

z ₀ is the position of the complex source and also determines the beam waist [24]. κ is the wavenumber; κz ₀ = 3.0 has been used throughout, corresponding to a beam with waist radius λ/2.

P_{n}^{m} (cos θ)

is an associated Legendre polynomial and γ_n,m is given by:

γ_{n, m} = {[\frac{(2 n + 1) (n - m)!}{4 π n (n + 1) (n + m)!}]}^{\frac{1}{2}}

2.2. Discrete Dipole Approximation

The DDA method has been reviewed elsewhere [25, 26], and details of our implementation were published previously [9]. To summarize, the scattering object is represented by an array of polarizable point dipoles, referred to as cells, situated on a cubic lattice of spacing δ. The field acting on any particular cell is given by the sum of the incident field and that scattered from each of the other cells. The field scattered by each cell is dipolar, with a strength determined by the product of the total field acting on the cell and its polarizability. These relationships may be written as a large, dense set of linear equations that are solved iteratively to give the polarization of each cell.

Traditional accounts of the DDA assume the scattering bodies are in vacuo. Optical trapping experiments, on the other hand, are usually performed in a medium such as water. By the scale invariance rule, the optical scattering for a fixed geometry depends on the wavenumbers in the different materials rather than the absolute values of the refractive indices [13]. Therefore the polarizability of the individual cells, α, may be written as [26]:

α = \frac{α_{CM}}{(1 - i \frac{2}{3} k^{3} α_{CM})} where α_{CM} = \frac{a^{3} (ɛ_{s} - ɛ_{m})}{(ɛ_{s} + 2 ɛ_{m})}

α has been corrected to ensure compatibility with the optical scattering theorem, and α _CM is the Claussius-Mossotti polarizability of the scatterer (permittivity ɛ _s), relative to the ambient medium (permittivity ɛ _m). a is the radius of a sphere whose volume is equal to one cell of the DDA lattice. This formulation may be extended to deal with optical anisotropy [27], in which case the polarizability of each DDA cell is represented by a symmetric tensor whose eigenvalues relate to those of the permittivity tensor through Eq. (5).

In the DDA, it is assumed that the optical force and torque on the particle may be evaluated by summing the contributions from each cell [6, 7]:

〈 F_{i}^{cell} 〉 = \frac{1}{2} Re (p_{j} \partial_{i} E_{j}^{★})

〈 T^{cell} 〉 = r \times 〈 F^{cell} 〉 + \frac{1}{2} Re (p \times E^{★})

Here p is the polarization of the cell obtained by solving the DDA equations and E is the total electric field at the cell. In Eq. (6), subscripts denote Cartesian components and repeated indices are summed. The angle brackets indicate cycle averages. The cell torque 〈T ^cell〉 consists of two contributions, an orbital part associated with the force on the cell, r × 〈F ^cell〉, and a spin part,

\frac{1}{2} Re (p \times E^{★})

, that arises due to optical anisotropy or absorption. Since the force is evaluated for a particle immersed in a medium, the values of p_i obtained from the DDA calculation, using the α of Eq. (5), must be multiplied by ɛ _m [23]. An alternative formulation of radiation forces on individual dipoles in the DDA is given by Hoekstra et al. [28].

As can be seen from this description, the method by which forces and torques are computed in the DDA is substantially different from that used in T-matrix theory (Section 2.1). Integration of an electromagnetic stress tensor, as is done in conjunction with the T-matrix method, is formally correct and can be found in numerous established text books (e.g. [23]). The method used in conjunction with the DDA, however, relies on the assumption that contributions from individual cells sum to give the required body forces. In fact, similar assumptions underpin the DDA approach per se, but convergence on continuum results has been rigorously justified [26]. That these two approaches are consistent is suggested in [29], for a slightly different set of equations. One of the purposes of the present study is to demonstrate the excellent level of agreement that may be obtained between these two, rather different, techniques, within the range of parameters encountered in optical trapping experiments. Furthermore, the method of calculating torques in the DDA, was published very recently [7]. To the best of our knowledge there are no other comparisons between torques calculated using this method and T-matrix calculations, either for non-spherical particles or for birefringent particles.

2.3. Stiffness coefficients

In general, a particle in a laser beam will experience forces and torques that are functions of its position and orientation. If a mechanical equilibrium configuration exists, the forces and torques may be linearly approximated in its vicinity [30]:

F = - K (q - q_{eqm}),

K = - {[\nabla F (q - q_{eqm})]}^{T} .

Here, F is the generalized force, (F_x, F_y, F_z, τ_x, τ_y, τ_z), and q the generalized coordinates (x,y,z, θ_x, θ_y, θ_z), where the vector θ = (θ_x, θ_y, θ_z) indicates a rotation of |θ| radians about an axis $\hat{θ}$ . q _eqm is the equilibrium trapping point. The stiffness matrix, K, decomposes into a 2 × 2 array of tensors:

K = [\begin{matrix} K^{tt} & K^{tr} \\ K^{rt} & K^{rr} \end{matrix}] .

in which the tensors K^tt and K^rr relate forces to displacements and torques to rotations, while the pseudo-tensors K^tr and K^rt couple rotations to forces, and displacements to torques, e.g.:

K_{x x}^{t t} = - \frac{\partial F_{x}}{\partial x}; K_{y y}^{r r} = - \frac{\partial τ_{y}}{\partial θ_{y}}; K_{x y}^{t r} = - \frac{\partial F_{x}}{\partial θ_{y}}; K_{y x}^{r t} = - \frac{\partial τ_{y}}{\partial x}

As previously discussed [30] the forces are generated by an optical flow rather than a potential. Consequently K need not be symmetric. This has implications for the thermal motion of the trapped particle, and for force and torque calibration. In general, K describes the behavior of forces close to equilibrium and, together with the hydrodynamic resistance, determines the stability of the equilibrium configuration. It plays an essential role whenever optical traps are used quantitatively. Because of its importance, and because it captures the local characteristics of the force field, we use it as a basis for comparing T-matrix and DDA calculations.

2.4. Qualitative description of optical forces

There are a variety of ways of describing optical forces. In the small particle limit, forces decompose into scattering and gradient components [31]. Larger, weakly scattering particles can be seen to behave in a quasi-dielectrophoretic manner, in that they seek to maximize their overlap with the incident intensity distribution. This behavior is modified by the presence of radiation pressure which, in general, pushes particles downstream but can also modify optical torques in non-trivial ways.

As noted above, the forces and torques on trapped particles may be evaluated from appropriate surface integrals of the Maxwell stress tensor. Using the divergence theorem, these may be converted to volume integrals, whose integrands can be interpreted as force and torque densities [22, 32]. For non-magnetic materials, these are given by:

〈 f 〉 = \frac{1}{2} Re [(J \times B^{*}) - \frac{1}{2} E_{j} E_{k}^{*} \nabla ɛ_{j k}],

〈 t 〉 = r \times 〈 f 〉 + \frac{1}{2} Re (D \times E^{*}) .

The force density 〈f〉 consists of two terms. The first relies on the electric current, J, and therefore depends on optical absorption to operate. The second includes the gradient of the permittivity, ∇ɛ_jk, and incorporates the effects of inhomogeneity and material boundaries. For an optically isotropic, non-absorbing dielectric, the first term vanishes and the force density reduces to:

〈 f 〉 = \frac{1}{4} {| E |}^{2} \nabla ɛ δ (n) \hat{n}

i.e. the force density is concentrated on the surface and is normal to it.

The torque density 〈t〉 also contains two terms. The first, r × 〈f〉, depends on the distribution of 〈f〉 over the particle. For a non-absorbing, optically isotropic material the fact that the force is normal to the surface implies that the effective body torque must be zero parallel to an axis of cylindrical symmetry i.e. an optical field cannot apply a torque to a sphere of this sort, while an optically isotropic cylinder diffuses freely about its axis. The second term, 〈D × E〉, depends on optical anisotropy. For optically isotropic media, it will be zero. Otherwise there will be a tendency to align the greatest component of the refractive index tensor with E.

The force and torque densities, 〈f〉 and 〈t〉, should not be confused with the contributions arising from individual cells in the DDA simulation, 〈F ^cell〉 and 〈T ^cell〉 [Eqs. (6) and (7)]. Although both sets of expressions may be integrated to yield the total body force and torque, they clearly represent different quantities and have different properties. For example, by Eq. (14), 〈f〉 is confined to the surface of an optically isotropic, non-absorbing dielectric and is normally directed. On the other hand, 〈F ^cell〉 need not be zero for a DDA cell in the interior of a modelled object nor, for a boundary cell, need it be normal to the surface on which it lies. It is one of the aims of the present study to demonstrate that the sums of 〈F ^cell〉 and 〈T ^cell〉 over the DDA cells do indeed converge on the theoretical values of force and torque as given rigorously by the T-matrix theory, and that this convergence can be reached without incurring excessive computational cost, for arbitrary particles of colloidal size and medium relative refractive index.

3. Numerical Results

In the following calculations, we use CSP beams for both the T-matrix and DDA results. The coordinate frame is shown in Fig. 1. Right handed coordinates are used with the origin on the focal point, the z-axis coincident with the beam axis, and the x-axis parallel to the polarization direction. Propagation is in the positive z-direction and the ambient medium is water, with refractive index n _m = 1.333, in which the laser wavelength is 0.6μm.

Fig. 1 Schematic showing the coordinate axes, and (a) a trapped sphere, (b) a prolate spheroid and (c) a Janus sphere.

Download Full Size | PDF

3.1. Dielectric spheres

For a dielectric sphere, the stiffness matrix is purely translational i.e. only K^tt is non-zero. This tensor is diagonal in the indicated coordinate frame (Fig. 1), and contains three independent stiffness coefficients, k_x, k_y and k_z, relating to restoring forces parallel to each of the axes. Figure 2 shows the results of DDA calculations of the trapping height, z _eqm, and the stiffness coefficient, k_x, for dielectric spheres of radii, r = 100, 200,300 and 400nm and refractive index n _s = 1.1n _m, as a function of the DDA lattice parameter, δ. The results for large δ are erratic, but ultimately settle on particular values of z _eqm and k_x when δ ≲ 25 nm. The most substantial error occurs in z _eqm for the smaller spheres. Based on these results, we use δ = 15nm in the work presented below, unless stated otherwise.

Fig. 2 (a) The equilibrium trapping height, z _eqm, and (b) stiffness, k_x, for dielectric spheres of radii r = 100,200,300 and 400nm, as a function of the DDA lattice parameter, δ. z _eqm is measured from the focal plane.

Download Full Size | PDF

The number of dipoles, N, used in Fig. 2 scales as expected with r ³ and with δ ⁻³. N ranges from 19, for r = 100 nm and δ = 60 nm to 267,731 for r = 400 nm and δ = 10 nm. The calculation times increase with N, the slowest being 64 s per force evaluation, when running in parallel on an 8 core, 2.5 GHz Xeon system. These times should be compared with a few seconds for a T-matrix calculation, or several hours for an equivalent FDTD simulation [5].

Figure 3 shows the effect of sphere radius on the equilibrium trapping height, z _eqm, and translational stiffness parallel to the polarisation direction, k_x. In Figs. 3a and 3b, the effect of lattice spacing, δ, is assessed. It can be seen that the overall shape of the curves is preserved for quite large values of δ, but that they become increasingly noisy for δ > 20nm. Figures 3c and 3d show the same curves for DDA calcuations with δ = 15nm compared with exact T-matrix calculations of the same system. The DDA and T-matrix results are so similar that they are visually indistinguishable. The greatest difference observed across the range of radii is c. 1%, while the average difference is c. 0.2%. Ambient media with three different values of n _m have been used confirming the need to incorporate ɛ _m in the DDA force calculation as described above (Section 2.2).

Fig. 3 (a) The equilibrium trapping height, z _eqm, and (b) the translational stiffness k_x, for dielectric spheres as a function of radius, for four values of the DDA lattice parameter, δ. The curves for δ ≥ 30nm are offset vertically, as they would otherwise be superimposed. (c) z _eqm and (d) k_x for dielectric spheres, comparing the DDA with δ = 15nm (black curves) and the T-matrix method (red curves), for three values of n _m.

Download Full Size | PDF

The results in Fig. 3 follow well recognized trends [33]. Geometric resonances cause the trapping height to oscillate as the sphere radius increases. The oscillation period corresponds, approximately, to a quarter of the wavelength of the incident light (or half a wavelength, in terms of the sphere diameter). When the sphere is very small, the momentum of the incident beam is only weakly deviated and the stiffness coefficients are small. On the other hand, when the sphere is large enough to encompass the focal region, the trap becomes less sensitive to the precise position of the particle and the stiffness coefficients are again low. Between these regimes the stiffness coefficients achieve maximal values, which occurs when the radius is approximately equal to half the wavelength in the ambient medium. This corresponds, in the present case, to the beam waist radius.

3.2. Geometric Anisotropy

When the dielectric spheres of the preceding sub-section are replaced with spheroids, the optical force field acquires a rotational aspect. Optical forces and torques will generally be present, whenever they are not prohibited on symmetry grounds [9]. In fact, optical torques can be applied to spheroids about any axis, apart from the symmetry axis [22]. For example, a prolate spheroid will be trapped beyond the focal point, in an upright orientation, with its long axis parallel to the beam axis [34]. Displacements parallel to the coordinate axes will be met by restoring forces. Similarly, rotations about either the x or y axes, but not the z axis, will give rise to restoring torques.

However, this is not the full set of optical forces and torques that act on the spheroid close to its equilibrium configuration. Displacements in the x or y directions also induce excess forces parallel to the beam axis, as well as coupling torques. The forces induced are even functions of displacement, while the coupling torques are odd functions. Similarly, rotations about the x and y axes produce excess axial forces (even in displacement) and coupling forces (odd in displacement). Forces that are even functions of displacement do not appear in the stiffness matrix. Therefore, many of the components of K vanish identically. The complete set of forces and torques acting on a prolate spheroid are analogous to those acting on a cylinder [30]. K^tt and K^rr will be diagonal, while K^rt and K^tr include coefficients that relate displacements in x and y directions to torques about the y and x axes, and rotations about the x and y axes to forces in the y and x directions, respectively.

In Fig. 4, we show the non-zero components of K, for prolate spheroids (n _s = 1.45) in water. The spheroids have an aspect ratio b/a = 2 and the results are plotted as a function of equivalent radius $r_{equiv}^{3} = a^{2} b$ . The lattice spacing used in the DDA calculations is δ = 10nm. Once again the DDA results are practically indistinguishable from the T-matrix results. The translational stiffness coefficients (Fig. 4a) behave similarly to those of a sphere: when the particle is too small to scatter strongly or big enough to envelope the focal region completely, the translational stiffness is low. Between these regimes, a maximum is reached. For the range of sizes considered, the rotational stiffness coefficients (Fig. 4b) increase monotonically with the size of the particle. This is due, in part, to the fact that rotation through a small angle, δθ, involves increasingly large displacements of the tip of the object, as the spheroid increases in height. This effect counterbalances the reduction in rotational stiffness one might expect to occur when the particle encompasses the focal region.

Fig. 4 Stiffness coefficients for prolate silica spheroids with an aspect ratio of two, as a function of equivalent radius. Calculations are performed using both the DDA and the T-matrix method.

Download Full Size | PDF

The coupling coefficients, Figs. 4c and 4d, derive from two sources: deviation in the linear momentum carried by the incident beam and the fact that trapping positions are displaced from the focal point. The first of these induces coupling when, for example, asymmetric reflection from an angled object gives rise to a force as well as a torque. The second arises on purely energetic grounds. Like spheres, spheroids are trapped beyond the focal point of the beam. When a spheroid is displaced in a transverse direction it will seek to increase its overlap with the incident beam profile both by reversing the displacement, and by rotating back towards the beam. Both of these effects tend to increase with size and, because of the role played by radiation pressure, they are necessarily non-conservative i.e. K is non-symmetric [30]. This is confirmed by Figs. 4c and 4d, from which it is clear that $K_{y x}^{r t} \neq K_{y x}^{t r}$ and $K_{x y}^{r t} \neq K_{x y}^{t r}$ .

Having established the accuracy of the DDA for geometrically anisotropic particles, more challenging systems may be investigated with confidence, i.e. shapes that would be difficult to simulate using a standard T-matrix approach, or very time consuming by FDTD. In Fig. 5, we examine the stability of colloidal silica particles in the shape of letters. The dimensions of the letters are shown in the figure. The optical trapping properties of an entire alphabet of such letters has previously been investigated experimentally [17]. Since the letters are flat, they tend to align in the plane of polarization (xz-plane). Therefore, attention is restricted to translation and rotation in this plane. The forces and torques experienced by the letters, in the vicinity of the experimentally established trapping configurations, are examined.

Fig. 5 The forces and torques acting on silica letters, computed using the DDA. The dimensions are: ‘N’, 4 × 5μm; ‘A’, 3 × 4.6μm; ‘O’, 3 × 5μm. All letters have a thickness of 0.5μm in the y-direction, and lie in the xz-plane. (a) The y component of torque, τ_y, when rotating a letter ‘N’ about the y-axis through the two centres shown. (b) τ_y for rotation of a letter ‘A’ about the y-axis through the centre of the letter. (c) The x component of force, F_x, for a letter ‘O’ being translated parallel to the x-axis, for the two orientations shown. (d) A schematic illustration of the preferred trapping orientations of the three letters.

Download Full Size | PDF

In Fig. 5a, the letter ‘N’ is rotated about an axis parallel to y through two pivot points, and the induced torque is plotted as a function of rotation angle, θ_y. The first pivot is located at the centre of the letter, on the diagonal, while the second is half way up the left hand vertical strut. As can be seen, rotation about the central pivot has a stable orientation at θ_y ≈ 0.2π rads. This corresponds to the diagonal strut of the ‘N’ being aligned with the beam axis. Rotation about the left hand pivot produces a stable orientation just off vertical (θ_y ≈ −0.05π rads). The two trapping orientations are illustrated schematically in Fig. 5d.

Figure 5b shows the induced torque on the letter ‘A’ as it is rotated about its centre. The legs of the ‘A’ are straight, as in [17]. In this geometry, rotational stability is achieved when the letter is lying on its side, i.e. θ_y ≈ ±π/2 rads (see Fig. 5d), as found experimentally.

Finally, Fig. 5c shows the x-component of the optical force as the letter ‘O’ is displaced in the x-direction. The solid curve is generated when the letter is oriented vertically, the dashed curve when its long axis points in the x-direction. As was observed experimentally [17], stable trapping only occurs in the horizontal orientation i.e. the stable configuration is with the beam passing through one of the smaller sides of the ‘O’, as illustrated in Fig. 5d.

3.3. Optical Anisotropy

Optical torque can also be generated by optical anisotropy. As in the case of a dielectric sphere, ∇ɛ will be zero inside a homogeneous birefringent sphere, and 〈f〉 will be confined to the surface and normal to it. As a result the boundary contribution to the torque vanishes identically. Therefore, 〈t〉 is entirely given by the second, spin-like, term in Eq. (13): 〈D × E〉. Since positive birefringence leads to the alignment of the optic axis with the polarization direction, this property allows control over the rotation of optically trapped birefringent particles [35].

In Fig. 6a, we plot the azimuthal component of the stiffness, $K_{z z}^{r r}$ , for rotations of a series of positively birefringent spheres about the beam axis. Since the variation of torque with respect to azimuthal angle is proportional to sin(2θ_z) [36], this stiffness is numerically equal to twice the available torque for rotation about the beam axis. Once again the DDA results reproduce the T-matrix calculations very closely. In general, the azimuthal stiffness increases both with birefringence and radius. However, 〈t〉 depends on the phase difference between E and D. For large path lengths the sign of 〈t〉 changes and so the magnitude of the torque, obtained by integrating 〈t〉 over the particle [22], diminishes, as can be seen in the figure for Δɛ = 0.4.

Fig. 6 (a) The azimuthal stiffness, $K_{z z}^{r r}$ , of homogeneous birefringent spheres in Gaussian beams as a function of sphere radius for several dielectric anisotropies, Δɛ. Calculations are performed using both T-matrix and DDA methods. (b) The azimuthal stiffness of quartz cylinders of various radii, as a function of length. The cylinders are oriented with long axes parallel to the beam axis. n _o = 1.544 parallel to the cylinder axis and n _e = 1.553 parallel to a radius. Calculations were performed using the DDA method only.

Download Full Size | PDF

The flexibility of the DDA allows other optically anisotropic structures to be modelled. In Fig. 6b we consider the angular behavior of optically trapped quartz cylinders. As with other elongated objects, the rods are trapped vertically, with their long axes aligned parallel to the beam. We consider the case where the optic axis is perpendicular to the cylinder axis, and is therefore preferentially aligned with the polarization direction, x. Such particles have applications in novel torque microscopes [18, 19]. In Fig. 6b we show the angular stiffness about the beam axis as a function of length, for quartz cylinders with several different radii and Δn = 0.009. The aspect ratios of the cylinders vary between c. 1.4 and 50. While increasing the length of the rod adds stability and increases the stiffness with which the cylinder is held vertically, it does not greatly improve the angular behavior about its symmetry axis. Indeed, it can be seen that there is a distinct change in behaviour between the narrower cylinders (0.1 ≤ r ≤ 0.2μm) and the wider ones (0.25 ≤ r ≤ 0.35μm). For the narrower cylinders $K_{z z}^{r r}$ , and hence the available optical torque, gradually increases with cylinder length. However, the rotational friction will increase approximately linearly with the length of the rod [37], thereby placing a limit on the maximum possible rotation rate at a given power. For the wider cylinders, the torque reaches a maximum for lengths in the range 2 – 3μm, and then decreases with further increases in cylinder length.

The explanation for the differing behaviour appears to relate to whether propagation of light down the cylinder length is permitted or prohibited. When the radius is below the waveguide cut-off, most of the torque will come from rays entering the rod from the sides. However, as soon as the radius exceeds this cut-off, propagation will occur down the length of the cylinder. Then, phase changes will build up between E and D, so that 〈t〉 goes negative and starts reducing the total torque i.e. the reduction in torque due to aggregated phase differences is only apparent when propagation down the rod is allowed. Consequently, decreasing n _o and n _e, while maintaining the same Δn, will lead to an increase in the cut-off radius and the reduction in torque will occur for wider rods, as has been observed (unpublished results).

It is clear, therefore, that the optimal torque will be obtained from the wider quartz cylinders considered, for lengths between 2 and 3 μm. Additional torque may be gained by increasing the radius even further, but this will become ineffective once the cylinder radius exceeds the beam waist radius. Interestingly, the quartz cylinders used in experiments have radii in the range 0.13 to 0.3 μm, and lengths around 1 μm, suggesting that there is scope for improving the torque that can be delivered experimentally by this method [18, 19]. Deufel et al. quote an azimuthal stiffness of 0.0114 pN.μm.rad.⁻¹mW⁻¹ for a cylinder with radius 0.265 μm and length 1.1 μm [18]. This is smaller than the values predicted in Fig. 6b, by a factor of 4. The difference may be attributed to a number of factors including differences in the refractive indices and laser wavelength, irregularities in the shapes of the fabricated cylinders, differences in the beam waist and differences in where the power is measured. The latter two points are likely to make the largest contribution to the observed discrepancy.

3.4. Inhomogeneity

Finally, we examine the affects of refractive index inhomogeneity. To test the accuracy with which the DDA can represent this, we examine the forces on core-shell spheres. We consider particles with a low index core (PTFE, n = 1.35) and a high index shell (silica, n = 1.45). In Fig. 7, we plot the equilibrium trapping height, z _eqm, and the three components of translational stiffness, k_x, k_y and k_z determined at z _eqm, comparing the results from DDA calculations and T-matrix theory. Good agreement is observed between the two methods for z _eqm (see Fig. 7b), but some discrepancies are found for the stiffness results (Fig. 7c), which are particularly obvious when the shell gets thinner (see Fig. 7d). For thin shells, the discrete nature of the DDA lattice appears as discontinuities in the gradients of the stiffness curves.

Fig. 7 (a) Schematic of the PTFE/silica core-shell spheres. (b) Equilibrium trapping heights and (c) trap stiffnesses for the core-shell spheres as a function of r_a. (d) Expanded view of the rectangular region in part (c), plotted as a function of the shell thickness, r_b – r_a. In all cases, the shell radius, r_b, is held constant at 0.75μm; the results from T-matrix and DDA calculations are compared.

Download Full Size | PDF

The optical force field experienced by core-shell spheres is more complicated than for homogeneous spheres. For example, z _eqm shows an obvious variation with core size. When the core is very small, the particle behaves like a silica sphere; when it is very big, the particle is effectively a PTFE sphere. For a range of intermediate core sizes, the shell traps preferentially and, as it is confined to the edge of the particle, so the particle is displaced in the trap. In particular, for core radii between 0.35 and 0.55μm, trapping through the silica shell is favored to the extent that there are two distinct trapping heights (see Fig. 7b). In the first of these, the centre of the sphere is downstream of the focus (z _eqm > 0) and the trap is concentrated on the part of the shell nearest to the focus. In the second, the centre of the sphere is upstream the focus (z _eqm < 0), and the opposite region of the shell is trapped. Evaluation of the trap stiffnesses indicates that, when z _eqm > 0, the trap is stable for all displacements whereas, for z _eqm < 0, it is unstable to displacements in x or y. Therefore, the trap stiffnesses shown in Fig. 7c are only given for z _eqm > 0.

When the shell thickness is such that it influences the trapping position, there is an increase in the vertical stiffness relative to the transverse stiffness (Fig. 7c). This is because that part of the shell that is responsible for trapping is actually a spherical cap (shown shaded in Fig. 7a), which has a reduced spatial extent parallel to the beam axis (vertically) compared with the transverse directions. Thus, by adjusting the shell thickness it is possible control the relative trap stiffnesses in different directions. This leads to the attractive proposition of being able to tailor the properties of the trap for particular applications. For example, according to Fig. 7c, a core radius of 0.6 μm would produce a quasi-isotropic trap with k_x ≃ k_y ≃ k_z.

The influence of the shell is observed in other directions. For example, the preference for trapping through the shell leads to an approximately ellipsoidal region in the xz-section through the force field, and a circular region in the xy-section (shown in black in Fig. 8) within which the force tends to zero. While the only stable equilibrium configuration is that described above, it is clear that there is a tendency for the core-shell particle to move away from the beam axis and away from the focus, in all directions, because of the need to overlap the shell with the more intense parts of the beam.

Fig. 8 (a) Longitudinal (xz-plane, y = 0) and (b) transverse (xy-plane, z = 0) vector force fields for a PTFE/silica core-shell sphere with an outer radius r_b = 1μm and core radius r_a = 0.33μm, in a Gaussian optical trap. The force amplitudes are highlighted through the background colour map; a black background indicates zero force. In (a) the trapping position is indicated by a cross.

Download Full Size | PDF

Another topical example of structural inhomogeneity is offered by biphasic Janus spheres [20]. Janus spheres consist of two distinct materials combined into a single sphere (see Fig. 1c). Since the force density is normal to each interface, including the internal one, these particles can experience an optical torque, despite being spherical. As with the core-shell particles this leads to complexities in the force field that are not present for homogeneous, isotropic spheres. In Fig. 9a we plot the torque about the y-axis, τ_y, as a function of rotation angle about the same axis, θ_y, for a range of Janus spheres with different radii. θ_y = 0 corresponds to the interface normal being aligned with the z-axis. As can be seen, the torques experienced are complicated in structure. As the sphere radius increases the numbers of maxima and minima also increase and the stability of the θ_y = 0 configuration switches. In fact, for most radii, the sphere is marginally stable to small rotations around θ_y = 0. However, for larger rotations, the configuration becomes unstable for all radii. For radii of 0.6 and 0.7 μm, stable rotations become apparent around θ_y = π/2 radians. These are accompanied by transverse forces, F_x (Fig. 9c) that act to move the higher refractive index part of the Janus sphere into the beam. The behaviour for rotation about the x-axis is slightly different (Fig. 9b). The configuration at θ_x = 0 is unstable to x-rotations for all sphere radii. However, the accompanying lateral force, F_y (Fig. 9d) is two orders of magnitude smaller than the previous case. These issues will be discussed at greater length in a forthcoming article.

Fig. 9 The torques acting on biphasic Janus spheres, positioned at the focal point of a Gaussian beam, as a function of rotation angle about (a) the y-axis and (b) the x-axis. The accompanying lateral forces are shown for rotations about (c) the y-axis and (d) the x-axis. Calculations were performed using the DDA method, for spheres with radii in the range 0.3 to 0.7 μm. The refractive indices used were n ₁ = 1.05n _m and n ₂ = 1.1n _m.

Download Full Size | PDF

4. Discussion & Conclusions

Results have been presented for the trapping positions, trap stiffnesses and optically-induced forces and torques, for a range of dielectric particles with differing symmetries. Excellent agreement has been demonstrated between results obtained with the DDA and exact T-matrix theory for dielectric spheres, geometrically and optically anisotropic particles, and also inhomogeneous core-shell particles. To the best of our knowledge, these are the first comparisons of DDA and T-matrix calculations for optical forces and torques on the lower symmetry classes of particle discussed. The results confirm that the DDA approach is reliable for each of the classes of particle considered. However, the limitations of a lattice-based approach have also been highlighted, with minor deviations from the T-matrix results being noted in the case of core-shell particles, when the thickness of the shell approaches the period of the DDA lattice.

The DDA method has a number of advantages over the T-matrix theory. Probably the biggest of these is the flexibility of the lattice-based approach, which means that a wide range of complex shapes can be modelled easily, without requiring adaptation of the basic algorithm. This latter point has been illustrated by examining more general systems including letter-shaped colloidal particles, biphasic Janus spheres and birefringent cylinders. Although our treatment is far from comprehensive, several new phenomena have been exposed which admit to further analysis. The complex morphology of the optical force fields associated with the colloidal letters gives rise to equilibrium configurations that are not intuitive and often not unique for a particle. Similarly the transition in the behaviour of the quartz rods as their radius crosses the waveguide cut-off, and the coupled forces and torques generated by biphasic particles, was not foreseen. Each of these effects warrant closer inspection. Although it would be possible to apply T-matrix theory in these cases, some effort would be required to generate the T-matrix for each class of particles being examined. In addition, the reliance of T-matrix theory on central VSWF expansions of the optical field has implications for the degree of geometric anisotropy that can be represented accurately [15]. In the DDA, by contrast, the scattered field is expanded about a set of points distributed over the particle, and does not suffer from this problem

It has been suggested elsewhere that the DDA method is very expensive in terms of computational resources. For example, one recent study reported that a model containing c.10⁴ dipoles took several hours to solve for the dipole moment on an 8-core Linux workstation [11]. Although our own implementation of the DDA method is slower than T-matrix theory, it is at least two orders of magnitude faster than the equivalent FDTD calculations, and the computational overheads incurred are not excessive. For example, the most complex shapes considered here were the letter-shaped colloidal particles of Section 3.2. These letters were modelled with a few hundred thousand dipoles: ‘N’, for example, contained 329,365 dipoles with a lattice spacing of δ = 20nm. The typical time required for each force calculation was ≲ 120 s, on an 8 core, 2.5 GHz Xeon system. The calculation ran comfortably within the 16GB of memory available on the system. This performance compares very favourably with reports in the literature [11].

It is worth noting that around 20 separate force calculations are required to establish the position of z _eqm for each combination of parameters considered above, and a further five evaluations are then needed to determine each stiffness coefficient. Therefore, it is vital to ensure that the DDA calculations are performed efficiently. We perform the DDA matrix inversion iteratively using a Krylov method [9], with the evaluation of matrix-vector products accelerated through the use of FFTs. Indeed, the use of FFTs is crucial to obtaining a rapid solution.

With these points in mind, it is clear that the DDA provides a rapid method for computing the optically-induced forces and torques on irregularly shaped particles. As demonstrated above, lowering the particle symmetry, or introducing inhomogeneities, e.g. in the case of Janus particles, gives rise to complexities in the force and torque fields that influence trap stability and can indicate the existence of multiple equilibrium configurations as well as possible nonconservative effects. A deeper study of these systems will be the subject of a future publication.

Acknowledgments

The authors are grateful to Research Councils UK for financial support under the Basic Technology Programme: A Dynamic Holographic Assembler translation grant. This work was carried out using the computational facilities of the Advanced Computing Research Centre, Bristol University: http://www.bris.ac.uk/acrc/.

References and links

1. E. Fällman and O. Axner, “Influence of a glass-water interface on the on-axis trapping of micrometer-sized spherical objects by optical tweezers,” Appl. Opt. 42, 3915–3926 (2003). [CrossRef] [PubMed]

2. Y. Harada and T. Asakura, “Radiation forces on a dielectric sphere in the Rayleigh scattering regime,” Opt. Commun. 124, 529–541 (1996). [CrossRef]

3. D. A. White, “Numerical modeling of optical gradient traps using the vector finite element method,” J. Comp. Phys. 159, 13–37 (2000). [CrossRef]

4. R. C. Gauthier, “Computation of the optical trapping force using an FDTD based technique,” Optics Express13, 3707–3718 (2005). URL http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-10-3707. [CrossRef] [PubMed]

5. D. Benito, S. H. Simpson, and S. Hanna, “FDTD simulations of forces on particles during holographic assembly,” Opt. Express16, 2942–2957 (2008). URL http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-5-2942. [CrossRef] [PubMed]

6. P. C. Chaumet, A. Rahmani, A. Sentenac, and G. W. Bryant, “Efficient computation of optical forces with the coupled dipole method,” Phys. Rev. E 72, 046708 (2005). [CrossRef]

7. P. C. Chaumet and C. Billaudeau, “Coupled dipole method to compute optical torque: Application to a micro-propeller,” J. Appl. Phys. 101, 023106 (2007). [CrossRef]

8. D. Bonessi, K. Bonin, and T. Walker, “Optical forces on particles of arbitrary shape and size,” J. Opt. A: Pure Appl. Opt. 9(8), S228–S234 (2007). [CrossRef]

9. S. H. Simpson and S. Hanna, “Holographic optical trapping of microrods and nanowires,” J. Opt. Soc. Am. A 27, 1255–1264 (2010). [CrossRef]

10. L. Ling, F. Zhou, L. Huang, and Z.-Y. Li, “Optical forces on arbitrary shaped particles in optical tweezers,” J. Appl. Phys. 108(7), 073110 (2010). [CrossRef]

11. V. L. Y. Loke, M. P. Mengüç, and T. A. Nieminen, “Discrete dipole approximation with surface interaction: Computational toolbox for MATLAB,” J. Quant. Spectrosc. Radiat. Transf. 112, 1711–1725 (2011). [CrossRef]

12. T. A. Nieminen, H. Rubinsztein-Dunlop, N. R. Heckenberg, and A. I. Bishop, “Numerical modelling of optical trapping,” Comput. Phys. Commun. 142, 468–471 (2001). [CrossRef]

13. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption and Emission of Light by Small Particles (Cambridge University Press, 2002).

14. S. H. Simpson and S. Hanna, “Numerical calculation of inter-particle forces arising in association with holographic assembly,” J. Opt. Soc. Am. A 23, 1419–1431 (2006). [CrossRef]

15. F. M. Kahnert, “Numerical methods in electromagnetic scattering theory,” J. Quant. Spectrosc. Radiat. Transf. 79, 775–824 (2003). [CrossRef]

16. V. L. Y. Loke, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “T-matrix calculation via discrete dipole approximation, point matching and exploiting symmetry,” J. Quant. Spectrosc. Radiat. Transf. 110, 1460–1471 (2009). [CrossRef]

17. J. N. Wilking and T. G. Mason, “Multiple trapped states and angular Kramers hopping of complex dielectric shapes in a simple optical trap,” Europhys. Lett. 81, 58005 (2008). [CrossRef]

18. C. Deufel, S. Forth, C. R. Simmons, S. Dejgosha, and M. D. Wang, “Nanofabricated quartz cylinders for angular trapping: DNA supercoiling torque detection,” Nat. Methods 4, 223–225 (2007). [CrossRef] [PubMed]

19. B. Gutierrez-Medina, J. O. L. Andreasson, W. J. Greenleaf, A. Laporta, and S. M. Block, “An optical apparatus for rotation and trapping,” Methods Enzymol. 474, 377–404 (2010). [CrossRef]

20. I. Kretzschmar and J. H. K. Song, “Surface-anisotropic spherical colloids in geometric and field confinement,” Curr. Opin. Colloid Interface Sci. 16, 84–95 (2011). [CrossRef]

21. Z. F. Lin and S. T. Chui, “Electromagnetic scattering by optically anisotropic magnetic particle,” Phys. Rev. E 69, 056614 (2004). [CrossRef]

22. S. H. Simpson and S. Hanna, “Optical angular momentum transfer by Laguerre-Gaussian beams,” J. Opt. Soc. Am. A 26, 625–638 (2009). [CrossRef]

23. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, 1941).

24. C. J. R. Sheppard and S. Saghafi, “Electromagnetic Gaussian beams beyond the paraxial approximation,” J. Opt. Soc. Am. A 16, 1381–1386 (1999). [CrossRef]

25. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]

26. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106, 558–589 (2007). [CrossRef]

27. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]

28. A. G. Hoekstra, M. Frijlink, L. B. F. M. Waters, and P. M. A. Sloot, “Radiation forces in the discrete dipole approximation,” J. Opt. Soc. Am. A 18, 1944–1953 (2001). [CrossRef]

29. B. T. Draine and J. C. Weingartner, “Radiative torques on interstellar grains.1. Superthermal spin-up,” Astrophys. J. 470, 551–565 (1996). [CrossRef]

30. S. H. Simpson and S. Hanna, “First-order nonconservative motion of optically trapped nonspherical particles,” Phys. Rev. E 82, 031,141 (2010).

31. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, “Observation of a single-beam gradient force optical trap for dielectric particles,” Opt. Lett. 11, 288–290 (1986). [CrossRef] [PubMed]

32. I. Brevik, “Experiments in phenomenological electrodynamics and the electromagnetic energy-momentum tensor,” Phys. Rep., Phys. Lett. 52, 133–201 (1979).

33. A. B. Stilgoe, T. A. Nieminen, G. Knoner, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “The effect of Mie resonances on trapping in optical tweezers,” Opt. Express16, 15039–15051 (2008). URL http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-4-2661. [CrossRef] [PubMed]

34. S. H. Simpson and S. Hanna, “Optical trapping of spheroidal particles in Gaussian beams,” J. Opt. Soc. Am. A 24, 430–443 (2007). [CrossRef]

35. M. E. J. Friese, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Optical alignment and spinning of laser-trapped microscopic particles,” Nature 394, 348–350 (1998). [CrossRef]

36. S. H. Simpson, D. C. Benito, and S. Hanna, “Polarization-induced torque in optical traps,” Phys. Rev. A 76, 043408 (2007). [CrossRef]

37. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).

Application of the discrete dipole approximation to optical trapping calculations of inhomogeneous and anisotropic particles

Abstract

1. Introduction

2. Methods

2.1. T-matrix theory

2.2. Discrete Dipole Approximation

2.3. Stiffness coefficients

2.4. Qualitative description of optical forces

3. Numerical Results

3.1. Dielectric spheres

3.2. Geometric Anisotropy

3.3. Optical Anisotropy

3.4. Inhomogeneity

4. Discussion & Conclusions

Acknowledgments

References and links

Cited By

Figures (9)

Equations (14)

Optics Express