1. Introduction

The first use of optical forces for trapping micron sized particles was presented by Ashkin in 1970 [1], who demonstrated the manipulation of latex spheres in water. From this seminal paper many diverse applications for optical trapping have been forthcoming from both biological and physical sciences. Holographic optical tweezers are a more recent innovation that use computer generated holograms to produce many concurrent optical traps [2, 3]. Holographic systems with more than 170 traps have been reported [4]. Holographic systems offer an advantage over more traditional scanning systems, in that they allow the trapping position of a particle to be controlled both parallel and perpendicular to the beam axis direction. By using several traps simultaneously, it is possible to use holographic tweezers as a micro- or nanoassembler. The potential applications of such a dynamic holographic assembler (DHA) are numerous, including for example, the “bottom-up” construction of nanomachines and photonic structures.

A variety of different beams may be used to bring particles into contact in the DHA, including linearly and circularly polarized Gaussian light fields [1] and optical vortices [5, 6]. The exact choice of beam will be determined by the need to translate or rotate the particles, which might be of arbitrary shape. The success of the DHA will depend on the degree of control which is possible over the forces and torques in the system.

In this paper the application of linearly and circularly polarized optical traps in holographic assembly will be explored using computer simulations. There are three main features of these beams which give rise to forces and torques on trapped particles: radiation pressure, the gradient force and polarization effects. Our aim here will be to calculate the magnitude of the forces and torques on particles trapped in multiple focused laser beams, using numerical modeling techniques.

There are many numerical methods available for solving electrodynamic problems. For particles much larger than the wavelength of light a geometric approximation can be taken whereas, for relatively small particles, the Rayleigh approximation can be used. Where the particle size is similar to the wavelength of light only direct solutions of Maxwell’s equations are suitable. Some methods such as finite element modeling [7], the separation of variables method [8] and the T-matrix method [9, 10, 11, 12, 13] have been used in this regard. Unfortunately, the finite element method scales poorly with simulation volume making it unsuitable for the present application. On the other hand, the separation of variables and T-Matrix methods can be relatively fast, but are best suited to high-symmetry systems.

In the present work, the finite-difference, time-domain (FDTD) method appears most suitable, as it scales linearly with simulation volume and can handle arbitrary structures. To date, FDTD has received only modest attention in the field of optical tweezers [14, 15, 16] because it requires considerable computing power to simulate a reasonably sized sample volume. As a result, previously published work has been limited either to small or to 2-dimensional models of isotropic dielectrics. However, the FDTD method is very powerful, with the potential to model multiple, arbitrarily shaped particles in a range of different traps. In this paper we discuss a full 3-dimensional FDTD approach to optical trapping, including a treatment of anisotropic particles, which runs on a multiprocessor supercomputer.We begin by outlining the method in some detail, and then proceed to apply it to a range of structures interacting with one or several Gaussian traps.

2. Method

The FDTD method was used to compute the electric (E) and magnetic (H) fields propagating in a system of dielectric particles and optical traps. In FDTD simulations the E and H fields are updated by numerical integration of Maxwell’s curl equations in the time domain. The simulation space is divided into a lattice of discrete cells which contain the E and H field components. In our implementation we employed the cubic Yee cell [17], in which the E components are positioned on the edges of the cell and the H components at the centers of the faces. The simulation progresses by leapfrogging the solutions for E and H in time and space. Different objects are represented on the grid by specifying the refractive indices of the appropriate Yee cells. Two types of non-absorbing dielectrics were used in the simulations: isotropic and negatively-birefringent, anisotropic materials. The background medium in each case was water. The isotropic material was given the refractive index of silica and the anisotropic material was modeled as quartz. The main parameters used in the simulations are summarized in Table 1. The wavelength chosen matches that of our holographic tweezer system. Unless stated otherwise, we consider the beam to propagate parallel to the z-axis of the simulation box, and to be polarized in the x-z plane.

2.1. Sourcing the electric and magnetic fields

The simulations in this report model beams with a Gaussian intensity distribution. The sourcing of this light field is critical to the correct operation of the optical trap as the forces on the particles depend on the distribution of the E and H fields in the simulation space. The conventional approach to sourcing an FDTD simulation is to update the E and H field components across a plane in such a way that the required field is reproduced on one side of it. The most frequent choice of source is a plane wave, in which case this strategy results in the so-called ‘soft source’ condition [18].

Table 1. The parameters used to define the beam and particle characteristics.

View Table | View all tables in this article

The formal reason that the soft source can be applied, lies in the fact that any electromagnetic field is completely determined by the values that the tangential components of E and H take on the surface of any smooth plane that partitions space. This is manifested, for example, in the field equivalence theorem [19] and the orthogonality conditions of plane wave spectra [19, 20]. In each of these cases the fields on the empty side of the boundary plane are thought of as being generated by fictitious complex currents lying in the surface, which are given by:

where n̂ is the surface normal, and the incident fields, E ^inc and H ^inc, representing the Gaussian beam are known. In the present case, the beam propagates parallel to the z axis, and is sourced across a horizontal plane near the top of the simulation box, i.e. n̂=-k. The currents may then be incorporated into Maxwell’s curl equations:

where Δ is the cell side length, ω is the frequency of the beam and θ (x,y,z) is the phase of the current source at a particular point which will vary across the source plane in such a way as to create the desired beam. Clearly M ^inc is not physical but is applied in the same way as J ^inc:

The sine function is preferred to the cosine in Eqs. 2 and 3, as this minimizes the discontinuity in the source field at the origin, which is the point of largest current magnitude. From here the time-stepping equations relating B with E and D with H are generated in the usual way [18].

Thus, to source a Gaussian beam the appropriate E ^inc and H ^inc are required. Although the paraxial wave equation describes the propagation of Gaussian beams, in fact, solutions to the paraxial wave equation are not appropriate for use as a source in the present case. This is because the E ^inc and H ^inc fields derived from the paraxial wave equation are not strictly solutions of Maxwell’s equations, with the result that they do not propagate correctly in the FDTD scheme. One way of ensuring that E ^inc and H ^inc are solutions of Maxwell’s equations is to use expansions in sets of eigenfunctions of the vector Helmholtz equation. This guarantees that, at whatever point the expansions are truncated, the resulting fields remain solutions of Maxwell’s equations. Truncation of the expansion also yields fields that can be evaluated quickly and are easy to manipulate. Here, an expansion in vector spherical wave functions (VSWFs) is used:

where a_nm and b_nm are complex expansion coefficients, RgM _nm and RgN _nm are regularised VSWFs [21] and k is the modulus of the wavevector in the background medium. Since RgM _nmand RgN _nm are curls of one another the expansion for H ^inc is immediately given:

In the following we are principally interested in sourcing Gaussian beams, although the method discussed so far is more generally applicable. In this case, we use coefficients for a fifth order Davis beam [22]. For a beam polarized parallel to the x axis, we use:

where:

and s=1/kw ₀ where w ₀ corresponds to the beam waist radius. For a right- or left-circularly polarized beam, only the coefficients with m=+1 or m=-1, respectively, are required.

Figure 1 illustrates the use of the expansion described above. A longitudinal section is shown through the intensity distribution of a beam derived from the fifth order Davis equation (Fig. 1(a)) and, for comparison, a beam derived using solutions of the paraxial wave equation (Fig. 1(b)). In the paraxial case, the intensity distribution of the field becomes flattened near the source and the focus is not in the centre of the simulation space. This distortion arises because the paraxial approximation assumes that the plane wave components making up the beam have propagation directions very close to the beam axis. In reality this cannot be the case, because the beam is highly focused. Use of the VSWF expansion greatly reduces these artefacts, as shown in Fig. 1(a).

Where multiple traps are required in a simulation, these may be generated by linear super-position. It is assumed in all the cases considered below, where multiple beams are used, that those beams are mutually coherent.

2.2. Calculation of forces and torques

The force and torque on the particle are calculated using the total field which comprises incident and scattered radiation produced by the interface of the particle and the background medium (water). The force and torque on the particle are given by:

 

Fig. 1. Side views of a Gaussian beam produced using (a) the VSWF expansion and (b) the paraxial approximation. The dashed line indicates the source plane and the shaded region indicates the position of a uniaxial perfectly matched layer (UPML), used to eliminate reflections from the simulation boundaries [18].

Download Full Size | PDF

where T͂ is the Maxwell stress tensor, defined by:

ε and μ are the relative permittivity and permeability of the background medium, angle brackets indicate cycle averaged quantities, dS is the unit normal to the surface, r is the vector from the centre of the volume to the surface, and E and H are the real parts of the electric and magnetic field amplitudes. Some authors have used an alternative and more ad hoc approach to calculating forces and torques, involving use of the Poynting vector [14, 16]. We prefer to eschew such methods because the FDTD approach yields all of the fields required for a totally rigourous calculation, and the use of a lattice makes the surface integrals in Eqs. (8) and (9) particularly straightforward.

The surface of integration chosen in Eqs. (8) and (9) was a cuboid, consisting of an integer number of Yee cells. Numerical integration was performed with one sample point per cell. Since E and H do not coincide in the Yee cell, the magnitudes of the field components at each sample point were obtained by interpolation of the eight nearest values. The forces and torques were evaluated under conditions of steady state, which occurred after between 14 and 100 optical time periods, depending on the size of the system. For the purpose of the present study, convergence was taken to be when the incremental change in force was less than 1% i.e. |ΔF|/|F|<0.01, when averaged over the previous eight cycles.

2.3. Anisotropic materials

Anisotropic particles were simulated using a method introduced by Pereda [23]. Here the permittivity is no longer a scalar value but is represented by a 3×3 permittivity tensor, ε. The electric displacement field, D, is updated using the standard time stepping equations and then E is calculated using:

where ε=κ ^-1. All three components of D are needed to calculate a single component of E. However, only one of the components of D is evaluated at the same location as the corresponding component of E; the other components are distributed around the Yee cell. Therefore, the remaining components needed in Eq. (11) are interpolated, as shown here for E_x(i, j, k):

i, j and k are the grid indices of the cell. Similar relations hold for the other components.

2.4. The FDTD code

The FDTD simulations used our existing in-house code [24, 25], modified as described above. Simulations were performed on the University of Bristol supercomputer, Blue Crystal, which utilises 2.6GHz dual core AMD Opteron processors, with 2 Gbytes of RAM per processor core and Infiniband high speed interconnections. The FDTD code is parallelized using domain decomposition; inter-processor communication uses the MPI message passing interface. Typically, simulations were performed using 16 processor cores. Example timings are shown in Table 2. It can be noted that the simulation time increases a little less than linearly with the number of cells used, because larger models require a smaller fraction of the simulation time to be spent on inter-processor communication.

Table 2. Sample timings for FDTD simulations of a 1µm diameter silica particle in a Gaussian trap, using various lattice sizes. Each simulation consists of 2800 iterations i.e. 14 optical time periods, running on 16 processor cores.

View Table | View all tables in this article

As well as depending on the number of cells, the total simulation time will depend on the physical dimension of the sample volume, and the contents of the simulation box. The physical dimension determines the time for a wavefront to cross the simulation box, and is easily taken into account. However, the biggest variation in simulation time arises from differences in the structures being simulated; these can often lead to multiple reflections which can take many optical periods to settle. As a consequence, the longest simulation times experienced in the present study were around two days.

3. Results

We present calculations of the forces and torques on a range of particles possessing different degrees of anisotropy, when placed in one or more linearly or circularly polarized Gaussian beams. This is preceded by an examination of the optimum simulation conditions, in terms of the FDTD cell dimensions and the size of source plane used. The results which follow are normalised for beams with a power of 10mW, except where indicated.

3.1. Grid and source optimisation

The first case considered is that of a single silica sphere trapped in a linearly polarized Gaussian beam. This case has been treated previously by a number of workers [14, 16]. Here, we consider those factors affecting the accuracy of our calculations, namely the size of the FDTD cell and the dimensions of the source plane.

It is generally accepted that the accuracy of FDTD calculations will improve on reducing the grid cell size, as the discretized Maxwell’s curl equations tend to the continuous case. Also, the fidelity with which an arbitrary shape can be represented on a cubic grid will improve at smaller grid cell sizes, and the effects of numerical dispersion will decrease [18]. However, the computational expense increases with the number of cells (see Table 2), and so a compromise must be sought. Figure 2(a) shows the force acting on a 1 µm diameter silica sphere placed along the axis of a linearly polarized Gaussian trap, for several different cell sizes, Δ. It can be seen that the forces are similar for each cell size and that, once the cell size is less than λ/40, there is no substantial change in the values of the forces. Thus there appears little advantage in reducing the cell dimensions to less than λ/40.

 

Fig. 2. (a) The force on a 1 µm diameter silica sphere lying on the axis of a linearly polarized Gaussian beam, simulated with different FDTD cell sizes. Positive distances indicate positions closer to the source. (b) The force on a 1µm diameter silica sphere positioned 300 nm below the trapping point in a linearly polarized Gaussian beam calculated for several source area sizes. The cell side dimension was kept at λ/40 nm throughout.

Download Full Size | PDF

The size of the source plane which produces the beam is also a determining factor of the accuracy of that beam. Ideally, to produce a faithful representation of the Davies beam the simulation should be sourced across an infinite plane. As the biggest contribution to the electromagnetic field is produced at the centre of the source plane it is possible to truncate the source to give a finite simulation space without significantly changing the beam. Figure 2(b) shows the variation of the force on a particle in a linearly polarized Gaussian trap produced by sources of different sizes. It is clear that the variation in force is small for all source sizes investigated and tends to a limit for large areas. The average force shown in the figure is 0.9 pN, with a variation of ± 2.5%. The choice of source size for the remainder of this paper (29.2µm², i.e. 360×360 cells or 9λ×9λ), was felt to be a reasonable compromise between accuracy and computational expense. The largest source area considered required a threefold increase in computer time, for only a 2.5% improvement in the accuracy of the force.

3.2. Single trap stiffness

Using the optimum cell size and source dimension indicated above we have calculated effective spring constants for small displacements of a 1µm diameter silica sphere from its equilibrium trapping position in a Gaussian trap.

The results shown in Table 3 are for particles displaced parallel to the coordinate axes, where z is parallel to the beam axis and the beam is polarized in the x-z plane. The trap stiffness is slightly larger parallel to the x-axis than parallel to the y-axis. This appears to be because the irradiance of a general linearly polarized Gaussian beam is not circularly symmetric, even though the intensity will be [22]. The stiffness parallel to the beam axis is more than a factor of six smaller than in the perpendicular directions which means the particle will have a greater range of movement parallel to the beam due to Brownian motion. The lower spring constant parallel to the z-axis results from the fact that the light intensity gradient is much smaller along the beam axis than in its cross section.

Table 3. The spring constants in three perpendicular directions for a 1µm silica sphere in water trapped by a 10mW linearly polarized laser beam with a Gaussian intensity profile.

View Table | View all tables in this article

3.3. Two spheres in parallel traps

When assembling nanostructures, trapped particles must be brought close together and this affects their trapping positions due to the scattered radiation produced by each particle. The magnitude of this effect was investigated by performing FDTD simulations of pairs of 1µm diameter silica particles in parallel traps. Initially the equilibrium trapping position for a single particle in a single trap was found. The particle was fixed in this position relative to the trap, and an identical second particle and trap were introduced into the simulation. The two traps were mutually coherent. The excess force on each particle was calculated as a function of sphere separation (see Fig. 3).

 

Fig. 3. The excess, optically induced forces (F_x, F_y and F_z) arising between two 1 µm diameter silica spheres held at their trapping positions in two parallel Gaussian traps. The sphere separation is defined as the minimum gap between the two spherical surfaces.

Download Full Size | PDF

The excess force in the separation direction changes as a damped sinusoid with period equal to the wavelength of the laser. The period is understandable because the scattered radiation from each particle will interfere with the beam incident on the other particle either constructively or destructively depending on the relative phase. The magnitude of the excess force decreases with increasing separation due to the reducing contribution of the other beam and the difference tends to zero as the separation increases. Similar behaviour has been predicted using the T-matrix method [11]. It is worth noting that the magnitude of this effect is very small compared to the trapping forces shown in Figs. 2(a) and 2(b), only rising to 0.1 pN when the spheres are almost touching. Interestingly, similar calculations on pairs of 0.5µm diameter spheres predict excess forces that are an order of magnitude larger than those for the 1µm sphere, and of opposite sign. This difference arises because the smaller particles are more efficient Rayleigh scatterers, which scatter a large amount of the incident light beam in a perpendicular direction to the beam axis. The implication of this finding is clear: although forces due to scattered light are generally negligible compared with trapping forces, they can become significant when the sphere separation is less than ~ 200nm and, potentially, could interfere with the assembly processes.

3.4. A cylinder trapped in multiple beams

An important building block for use with the DHA is a dumbbell, which can be used to create 3-dimensional structures and tools for nano-manipulation. A dumbbell may be constructed using two silica spheres attached to the ends of a cylinder. In order to create such a dumbbell it is important therefore to understand the mechanism of cylinder manipulation using the DHA (see Fig. 4).

 

Fig. 4. Side view showing the relative intensity of three Gaussian beams trapping a 3µm long silica rod with radius 45 nm. The dashed line indicates the source. The intensity variation is shown as false colours varying from red (most intense) to purple (least intense). The vertical scale is stretched compared with the horizontal scale.

Download Full Size | PDF

In a single beam a dielectric cylinder will preferentially orient with its length parallel to the beam axis, as this maximises the amount of high refractive index material in the region of the beam with the greatest light intensity. Therefore, the orientation of a non-absorbing rod cannot easily be controlled using a single beam. However, the rod may be trapped and manipulated perpendicular to the beam axis if several beams are used [26]. In the following simulations, between two and five beams were used to trap a silica cylinder of length 3 µm and radius 45nm (see Fig. 4). For each simulation the results are normalised to a total beam power of 10mW; thus it is the distribution of light intensity along the cylinder that is being compared.

Effective spring constants were calculated for displacements of the cylinder parallel to each coordinate axis, for varying numbers of traps (see Table 4). The order of magnitude difference between the spring constants in the z and y directions results from the fact that the intensity gradient in the beam direction is much less than that in the cross section of the beam. Comparing the spring constants in the x-direction, it can be seen that the system has a higher stiffness when the laser power is concentrated towards the ends of the cylinder, i.e. when only two beams are used, decreasing as the number of beams increases.

Table 4. Effective spring constants for displacements parallel to the coordinate axes defined in Fig. 4, for a silica cylinder trapped in multiple Gaussian beams.

View Table | View all tables in this article

 

Fig. 5. The restoring force on a cylinder trapped by three parallel traps, as it is moved parallel to the x-axis from the central position seen in Fig. 4. The cylinder is 3µm in length, and has a radius of 45 nm. The traps are 1.14µm apart.

Download Full Size | PDF

The magnitude of the spring constant in the x direction is similar to that in the z direction, because only light incident close to the ends of the cylinder contributes to the gradient force in the x direction. As a result, the rod is able to slide relatively easily along its long axis. For example, Fig. 5 shows the restoring force as the cylinder is displaced along the x axis in a system with three traps. At zero displacement, there is no restoring force, corresponding to the position shown in Fig. 4. The restoring force is a maximum at a displacement of ~ 400 nm, as one trap is now at the end of the cylinder. Clearly, to reduce the extent of this sliding the outermost traps should be positioned as close as possible to the ends of the cylinder.

The effects of rotations on the cylinder are shown in Figs. 6 and 7. In Fig. 6, the cylinder was rotated about the y axis i.e. the axis perpendicular to the beam and rod major axis, and the resulting torque on the particle is shown. For angles less than 45°, the beams apply a restoring torque in all cases, indicating that the horizontal configuration is stable. In the case of the two beam simulation the torque remains negative at larger angles, as there is little light intensity between the beams. This differs from the three beam case where the central beam causes the torque to change sign when the cylinder is at about 45°, stabilizing the vertical orientation when the rod is rotated by more than 45°. The simulations with four and five beams show a smaller magnitude of torque because the light is distributed along the cylinder and the intensity gradients are much reduced.

 

Fig. 6. Torque on a dielectric cylinder as it is rotated about the y-axis, in varying numbers of parallel traps.

Download Full Size | PDF

 

Fig. 7. (a) The torque on a dielectric cylinder held in several parallel traps, as it is rotated about the beam axis. (b) Diagram showing the view along the beam axis for a cylinder in three beams. (c) Similar diagram for five beams.

Download Full Size | PDF

The torque on the cylinder when rotated about the beam axis is shown in Fig. 7(a). In all cases the maximum restoring torque occurs at approximately 10°. Two types of behaviour are observed. For 2 and 3 beams, the restoring torque drops off quite rapidly with angle, as the rod moves out of the range of the outermost traps (illustrated in Fig. 7(b)). For 4 and 5 beams. however, the drop in restoring torque is more gradual, as the inner traps exert an influence to larger angles (Fig. 7(c)). Overall the maximum restoring torque diminishes with number of traps, if the total power available is kept constant.

3.5. Oblate silica particles trapped in a circularly polarized beam

In a recent paper oblate spheroids trapped in a circularly polarized Gaussian beam were predicted to rotate about the beam axis in different directions depending on their size [12]. The results were initially produced using the T-matrix method; here we demonstrate the same phenomenon using FDTD simulations. A range of different sizes of oblate spheroids with an aspect ratio, d, in the range 0.2–0.5 as well as a flat, coin-shaped particle, were positioned at the trapping point of the beam with their short axes perpendicular to the beam axis, and the torque about the beam axis was calculated (see Fig. 8). The thickness of the coin-shaped particle was 150 nm. As can be seen in Fig. 8, all of the particles show torque reversal over some size range, apart from the oblate spheroids with δ=0.2.

 

Fig. 8. The torque exerted by a circularly polarized Gaussian beam on oblate spheroids with aspect ratios, δ, in the range 0.2-0.5 and a flat coin shape with thickness 150 nm.

Download Full Size | PDF

There are three distinct regions of particle size which affect the direction and magnitude of the rotation of the particle. For relatively small particle sizes where the particle is completely contained by the high intensity region of the beam, the torque is small. This is expected, because the torque increases with the volume of particle in the beam. The direction of rotation in this region is with the rotating polarisation and the rate at which the angular momentum, carried by the incident beam, is transferred to the particle is positive implying a steady loss of angular momentum from the beam itself. The second region has particle dimensions which are similar to the beam waist. The rotation is in the opposite sense to the polarisation. This implies that the angular momentum is being steadily transferred to the beam as a result of its interaction with the particle. Because the values taken by the spin component of the angular momentum are confined to h̄ω per photon, it is presumed that the interaction is mediated by an orbital component as was suggested in [12]. Scattering of light from particles whose symmetry is less than that of a sphere, generally results in fields whose azimuthal dependence implies the presence of orbital electromagnetic angular momentum. Particles in the third region are only partially contained within the beam and have high magnitudes of torque because they have large volumes. These particles rotate in the same direction as the beam. As the size is increased further the direction of rotation does not appear to change but increases in size approximately in proportion to the volume of the particle within the beam.

3.6. Trapping of anisotropic spheres

The anisotropic material used in this report was calcite which is optically uniaxial, and has a reasonably large negative birefringence (see Table 1). This material is of interest because it orients with its optic axis perpendicular to the beam axis, unlike materials with positive birefringence which orient with unique axis parallel to the beam. Studies were performed using both linearly and circularly polarized traps.

Figure 9(a) shows the variation in the torque on a 1µm diameter calcite sphere with the angle, ϕ, between the optic axis and polarization direction in a linearly polarized Gaussian trap. The angle, ϕ, is measured in a plane perpendicular to the beam axis, and the torque is calculated about the beam axis. The results show a sin2ϕ dependence, which compares favourably with previous work on geometrically anisotropic particles that made used of the T-matrix formalism [13]. However it is difficult to apply the T-matrix method to optically anisotropic particles (see e.g. [27, 28]), which highlights the versatility of the FDTD approach.

 

Fig. 9. (a) The torque on a 1µm diameter sphere of calcite trapped in a 10mW linearly polarized beam. Both calculated results (points) and a sin2ϕ fit (solid line) are shown. (b) The torque about the beam axis acting on different sizes of calcite sphere in a circularly polarized Gaussian beam.

Download Full Size | PDF

In a second study the calcite sphere was placed in a circularly polarized Gaussian beam to see if any reversal in rotation was observed in optically anisotropic systems, by analogy with the effects reported in Fig. 8. The results are shown in Fig. 9(b), which shows the induced torque about the beam axis as a function of particle size. It is clear from the figure that torque reversal does not occur for calcite, although we cannot rule it out for materials with different birefringences. For particles with dimensions smaller than the beam waist, w ₀, i.e. up to ca. 450nm in the figure, the torque scales approximately with the volume of the particle. However, for particles with radius greater than w ₀ the lateral edges of the particle, which are outside the beam, contribute little to the torque and so the torque varies linearly with the radius of the particle.

4. Discussion

The FDTD approach is clearly a versatile model for determining the interaction of light with matter and has the potential to be used to calculate forces and torques on quite complex arrangements of particles and traps. However, it is important to recognise the limitations of the method. The accuracy of the FDTD simulations is limited by the ability of the technique to represent an arbitrary particle faithfully using cubic cells. This automatically places a lower limit on the size of particles that may be simulated. It is always possible to increase the number of Yee cells in the calculation, i.e. to use more, but smaller cells. However, in 3-dimensions, the number of cells used rapidly becomes prohibitive, even when using a supercomputer, so it is important to find a compromise between computation time and accuracy. For the simulations presented here, a cell dimension, Δ, of 15nm or λ/40 was selected as optimal, on the grounds that the computed optical forces on 1µm diameter silica spheres appear to converge at this lattice dimension. The surface roughness arising from the use of the lattice in this case is around 1/60 of the sphere diameter and may be regarded as negligible. However, in the simulations of silica cylinders presented in Section 3.4, the use of a 15 nm cell dimension is quite crude, and the cylindrical cross section is far from circular. For the purposes of the present study, this inaccuracy was deemed not to be important.

We use Eqs. (1) to (7) to source a beam which propagates on the FDTD lattice, and which is Gaussian in three dimensions to a high degree of accuracy. However, the concept of a Gaussian beam is based on the paraxial approximation to the Helmholtz equation. It is arguable that in the current situation this approximation no longer pertains: laboratory beams need not be particularly close to paraxiality and, as a consequence, our beams may differ from those used experimentally. Effects such as the torque reversal predicted in Section 3.5 appear to depend quite subtly on the overlap between the spheroidal particle and the beam waist, and therefore may differ qualitatively as the beam shape changes. The sourcing of more accurate beams that take into account optical defects such as astigmatism, is a topic for future study.

There are three mechanisms which act on the particles that have been modeled: radiation pressure, intensity gradient forces, and torques due to the polarization of the beam. The radiation pressure and gradient forces determine the vertical trapping position of the particles and also the interparticle excess forces in the parallel beam simulations. The excess forces calculated are, in most cases, very small compared with the restoring forces found for small displacements. For example, given the trap stiffnesses quoted in Table 1 for a 1µm diameter particle, an excess force of 0.02pN corresponds to a horizontal displacement of 1 nm or a vertical displacement of 7 nm from the equilibrium trapping position. However, at very small separations (≲ 200nm), the excess force increases dramatically, and this effect is compounded when smaller particles are considered. At close separations, the excess optical force, F_x, can be comparable with the Van der Waals attraction, and may interfere with attempts to assemble structures. For example, taking the Hamaker constant for silica spheres in water to be A ₁₃₁=0.8×10^-20 J, gives an attractive force of 0.033pN for 0.5µm diameter spheres, or 0.067pN for 1µm diameter spheres separated by 100 nm, compared with an optically induced 0.08pN repulsion for the 1µm spheres and 0.17pN attraction for the 0.5µm spheres. Although we have assumed the trap geometry to be vertical, no account has been taken of the gravitational force on the silica spheres. In fact the force on a 1µm diameter sphere in water will be 0.007 pN. This will lead to a vertical change in trapping position of ~ 2 nm, and may safely be ignored.

The calculations based on the silica cylinder, highlight the importance of intensity gradient forces in manipulating extended structures. The intensity gradients provide the necessary restoring forces both to orient the cylinder and to localise its position in space. The use of multiple traps stabilises the horizontal orientation of the cylinder where it would otherwise tend to align coaxially with a single beam. It is important to place traps close to the ends of the cylinder when manipulating it, to ensure that the greatest intensity gradient is operating to prevent the cylinder from sliding along its length.

The sin 2ϕ dependance of torque on setting angle (Fig. 9(a)) provides a satisfying demonstration that optical anisotropy leads to the same behaviour as geometrical anisotropy when particles are held in linearly polarized traps, though not when circularly polarized traps are used. It also serves to validate the operation of our anisotropic FDTD code. That the orientation of anisotropic particles can be controlled by the beam polarisation has been well known for many years, and provides another means of orientational control that may be exploited when assembling structures in the DHA. The occurrence of torque reversal for geometrically anisotropic particles in circularly polarized beams has now been predicted using both T matrix and FDTD formalisms. The effect is quite subtle, depending on the overlap of the particle with the perimeter of the trap, but is large compared with any numerical errors occurring during integration of the Maxwell stress tensor. We are currently fabricating oblate spheroids with which to test the prediction experimentally. The absence of torque reversal for optically anisotropic particles is intriguing and is the subject of further studies.

5. Conclusions

We have demonstrated the application of a full 3-dimensional FDTD code to the trapping of optically isotropic and anisotropic dielectric particles in multiple Gaussian beams. The results presented are intended to provide guidance for the assembly of micro- and nanostructures using the DHA. For 1µm diameter spheres, we have shown that a FDTD grid cell size of 15nm or λ/40 is adequate for the calculation of accurate forces and torques. Finer grids will be needed when dealing with smaller systems. A method has been shown for sourcing the simulation with a Gaussian beam, by means of a multipole expansion of a fifth order Davis beam. The dimensions of the source plane are important in determining the absolute values of the forces calculated. By way of a compromise between computational expense and accuracy, we have standardised on a source plane measuring 9λ×9λ or 360×360 cells.

Simulations of two spheres in parallel traps indicate the presence of an additional, optically-induced force between the spheres, whose magnitude oscillates with sphere separation with a spatial wavelength equal to that of the incident beams. Although the magnitude of these forces are generally small compared with trapping forces, they become significant when the particles approach to within 200 nm, when they rival the interparticle Van der Waals forces. Simulations of a silica cylinder suggest that the use of multiple traps can stabilize the orientation of the cylinder perpendicular to the beam axis, but that for optimal control, the light should be concentrated towards the ends of the particle.

The simulations of spheroidal particles in circularly polarized beams confirmed the occurrence of torque reversal over a limited size range, as predicted previously using the T matrix method. In fact torque reversal is predicted for all particles considered except for the most anisotropic (δ=0.2). Optically anisotropic particles behave in the same manner as geometrically anisotropic particles in linearly polarized beams, displaying the same sin 2ϕ torque dependence. However, optically anisotropic particles do not appear to display torque reversal in circularly polarized beams at any size scale so far simulated.