We introduce a method whereby the electromagnetic field that governs the force on a Rayleigh particle can be tailored such that the resultant force field conforms to a desired geometry. The electromagnetic field is expanded as a set of vector spherical wavefunctions (VSWFs) that describe the field over all space. Given the incident field, the resultant force on a given Rayleigh particle can be calculated throughout a volume of interest. We use an evolutionary algorithm (EA) to search the space of coefficients governing the VSWFs for those that produce the desired force field. We demonstrate how Maxwell’s equations will support an “optical tunnel” that guides particles to a trap location while at the same time preventing particles outside the tunnel from approaching the trap. This result is of interest because the field is impressed throughout the domain; that is to say, once the field is generated, no additional control is required to guide the particles.
© 2011 OSA
Optical trapping of particles using tightly focused laser beams was first demonstrated by Ashkin et al.  and has since found use in many fields [2–4]. Recently, significant attention has been focused on structured light beams and in particular their application to optical tweezers . A number of groups have shown that arrays of traps can be generated and dynamically manipulated in real-time to guide and arrange particles as desired [6–12] with a good recent review given in . Alternate trapping beams have also been considered, most notably, “donut beams” composed of Laguerre-Gaussian modes where a high-intensity ring surrounds a low-intensity trapping region , non-diffracting beams with self-regenerating propagation characteristics [15, 16], and “optical bottle” beams composed of a dark central region surrounded in all directions by high intensity regions . The means by which these beams and traps are generated is varied, and includes multi-beam interference , the generalized phase contrast method , diffractive optical elements (DOEs) , holograms , and time-sharing traps [22–25].
Within the trapping literature, the desire to create arbitrary intensity distributions within a volume of interest has been strong, with recent work seeking to design three-dimensional (3D) fields [26–29]. Here, we are interested in sculpting 3D fields as well but would like to design fields that force particles through predefined paths in a volume of interest. We consider the possibility of creating 3D “light wires” with arbitrary geometry and demonstrate here that Maxwell’s equations will support the creation of an optical tunnel that guides particles within a defined region to a central trap and repels all particles outside the tunnel region.
We demonstrate that such a field is possible by expanding a 3D monochromatic field with a basis of vector spherical wavefunctions. These functions are each complete solutions to the vector Helmholtz equation and involve no paraxial approximation. We employ an evolutionary algorithm to search the space of complex coefficients governing the VSWFs and find a set of coefficients that generates a field closely matching a predefined field in a least-squares sense. While it is straightforward to directly compare desired electric fields, we take an additional step and search the space of force fields that a particular Rayleigh particle would experience if present in the field. Matching force fields governing Rayleigh particle dynamics assumes the effect of the particle on the impressed field is negligible, but we intend in future work to extend the method to particles that are on the order of the incident beam’s wavelength. The use of VSWFs allows for this possibility, which we discuss in Section 5.
In Section 2 we present the VSWF framework with which we parameterize 3D fields. In Section 3 we show how our numerical methods allow faithful reproduction of analytical results for a Gaussian beam and explain how the target force distribution was generated. We show how DE is used to search for desired sets of coefficients in Section 4, discuss results and directions for future study in Section 5, and conclude in Section 6.
2. Field expansion
We assume a homogeneous, absorption-free, source-free and dielectric medium and a monochromatic, time-harmonic electric field given byEq. (2), thus, any arbitrary complex amplitude can be written as a linear combination of the VSWFs as follows: 30].
The singularity-free regular VSWFs are given by30–33].
Choosing Nmax as the upper limit of the summation in Eq. (3) we can determine E(r) at any point within the volume of interest given values for the complex expansion coefficients and (for notational convenience, we combine the expansion coefficients into a single complex vector, ). Here we are not concerned with analysis, that is, determining the expansion coefficients given a defined field E(r) as discussed in . Rather, we synthesize E(r) given desired particle trajectories in the volume of interest.
3. Force field calculation
Given the electric and magnetic fields, the gradient and scattering forces associated with a given Rayleigh particle can be calculated. Following the development in , the instantaneous gradient force is given byEqs. (3), (4), (12), (15), (16), and (18).
Proper selection of the prefactor A in Eq. (3) is required to appropriately account for beam power. There has been some disagreement in the literature about how incident power, Pi, should be calculated for beams expanded as VSWFs. Moine and Stout  were careful to ensure proper normalization of beam power for particles with radii on the order of beam wavelength and questioned the idea that Pi ∝ ||ã||2, where the double-bars represent the Euclidean norm. This proportionality was later shown by Simpson and Hanna  to be true for a wide variety of laboratory beams. By properly accounting for A in Eq. (21) of  we calculate the prefactor as
We demonstrate the accuracy of the calculated fields by comparing forces calculated using the above procedure with the analytical results for a Gaussian beam presented by Harada and Asakura . The VSWF coefficients required to represent the Gaussian beam are found by using the development presented in . First, the analytic expressions for the plane-wave coefficients are given as36]. A Davis-type model of a Gaussian beam  allows for nonparaxial corrections by taking s 2 power series approximations of the vector potential, where the dimensionless beam shape parameter is defined by
The analytical description of the scattering force as calculated in  is given byEqs. (17) and (18) of . A moderately-focused, x-polarized (êi = [1 0]) beam acting on polystyrene spheres in water was assumed in  with physical parameters np = 1.592, nm = 1.332, w 0 = 5 μm, Pi = 100 mW, and a vacuum wavelength λ 0 = 0.5145 μm. We use the same values and calculate the Gaussian fields using Nmax = 130. Analytical and computational results for the Gaussian beam are compared in Figure 1.
We begin building the target force field, F T (r), by generating a Gaussian beam for a particle with R = 50 nm. Fewer coefficients (Nmax = 12) are required to accurately represent a more focused beam with w 0 = λm where λm = λ 0/nm. We discretize the z = 0 plane using increments of 0.2λm in the x- and y-directions and calculate the force field corresponding to a particle with the same properties used above to yield F g(r)|z =0. This field is then modified by leaving the central core of the field unchanged and inverting the sign of forces in the plane outside the central core as shown in Figure 2A. The target field in the z = 0 plane is then copied at additional z-planes to build the full 3D optical tunnel as shown in Figure 2B. Care must be taken to ensure that the spacing between sampled z planes is small enough to accurately capture all relevant spatial structure of the beam. Here, a spacing of 0.1λm between planes ensures that the computed gradient in the z-direction is accurate.
4. Design procedure
Our goal is to search the space of VSWF coefficients governing Eq. (3) for those coefficients that yield a force field that most closely matches the desired F T (r) in a least-squares sense. We combine all the coefficients into a single vector by letting33].
Given the large number of parameters governing the optimization (31) it would be difficult to implement a derivative-based optimization procedure like gradient descent because the partial derivative of (31) with respect to each of the parameters in v would need to be estimated numerically. For any reasonable value of Nmax this would quickly become prohibitive. Furthermore, such an algorithm is highly likely to converge to a local optimum if the search space is composed of multiple optima.
Because we cannot claim a priori that the search space is convex, we make use of a global direct search routine to search the parameter space, in particular, a specific type of evolutionary algorithm called differential evolution (DE) . Briefly, evolutionary algorithms are a class of optimization procedure that employ a population of randomly generated solutions and principles inspired by Darwinian evolution to efficiently search a parameter space for optimal solutions. The fitness (if maximizing) or cost (if minimizing) of each solution in the population is measured against an objective function and those solutions that best meet the objective are selected to “breed” and pass their characteristics on to the next generation of solutions. Here, we actually convert the minimization given by (31) into a maximization problem by taking the inverse of the least-squared error, but either approach will work. Over many generations the population will converge to an optimum solution. The stochastic nature of the search and the presence of many solutions in the population helps the algorithm efficiently search the space while avoiding convergence to locally optimal solutions.
DE is designed to work with vector solutions which, in this case, are the parameter vectors, v, governing E(r). During each generation, the algorithm steps through each member of the population and compares the fitness of that member (the “target vector”) to the fitness of a new vector that is created from the population (the “trial vector”). Each trial vector is created in the following manner: first, a “difference vector” is created by taking the difference between two randomly selected members of the population and multiplying by a scale factor Fs. Next, a third member of the population is randomly selected and added to the scaled difference vector to yield a composite vector. The final trial vector is created by populating each element of the vector with an element from either the target vector or the composite vector where an element from the target vector is selected with probability CR. Subsequently, the fitness of the trial vector is calculated and compared to the fitness of the target vector with the higher-fitness vector selected as a member of the next generation population. In this work we use Fs = 0.9 and CR = 0.5 as recommended in , but the optimizations are not especially sensitive to the exact values of these parameters. Once the next generation population has been created, the process repeats until a stop criterion is reached. Here, the algorithm runs for a fixed number of generations (500).
5. Results and discussion
The DE algorithm is capable of producing a solution that satisfies the main requirements of the desired optical tunnel. In particular, a field was designed with a core tunnel region that prevents particles that begin outside the core region from entering the tunnel. A sample force field from the best discovered solution out of 10 separate optimization runs is shown in Figure 3. Figure 6 shows the three Cartesian components of the electric field for three sample z planes. All of the runs converged to similar solutions.
The discovered force field is characterized by alternating rings of attraction and repulsion. Within the central region, particles are pushed toward the central axis and at some radius, particles are no longer drawn toward the axis but are instead pushed away. At still larger distances from the central axis, the force again switches direction and pushes particles back toward the center. These alternating bands can be seen as faint rings within Figure 3A.
Forces were only matched at the displayed points. In order to determine whether the field behaves similarly for points that were not explicitly matched during the optimization, we track a series of trajectories throughout the domain. Several trajectories–some beginning within the core region, some outside–were initiated from the plane at z/λ = −1. Each particle was always assumed to instantaneously move at terminal velocity within the fluid (assuming Stokes drag) and a finely sampled cube of points centered at the particle’s current location was used to determine the local forces present on the particle. The results of this exercise are shown in Figures 4 and 5.
Figures 4 and 5 illustrate that the algorithm has converged to a solution that is tunnel-like over most of the observed volume but is actually composed of a trap within the core (near the z/λ = 1 plane) where particles outside the core region circle about the center. The algorithm seeks to maximize the objective function over the optimization region but can only do so in a least-squares sense under the constraints of Maxwell’s equations. A more tunnel-like field might be possible if the size of the particle were increased such that the scattering force contributed more to the overall field. In addition, a longer tunnel region could be generated if the target field were extended, but at the expense of increased computation time. Regardless, the present field does succeed in generating a core region that is protected from particles that originate beyond the core boundary.
Figure 6 shows the progression of the z-component of the Poynting vector at various points along the z-axis as well as all components of the Poynting vector at the z = 0 plane. The structure of the discovered beam is quite different from both Laguerre-Gauss  and Bessel beams . Note that the intensity profile of the beam changes drastically as it propagates along the z-axis. In particular, the relative strength of the high-intensity rings changes significantly, which is not typical of Bessel or Laguerre-Gauss beams. The direction of the Poynting vector also oscillates at points located on the z-axis, which is another differentiating feature of the discovered beam.
The main goal of this work, however, is to illustrate a procedure by which optical geometries can be designed. In particular, searching the space of VSWF coefficients may be an efficient means of determining whether a given field can be produced by a diffractive optical element (DOE). DOEs have been used for some time to sculpt light into desired patterns , however, it is generally not known a priori whether the desired field is consistent with Maxwell’s equations . The algorithms used to produce the DOE, such as Gerchberg-Saxton (GS) , will converge to some approximate solution, but it is not possible to determine whether there is a better approximation. In addition, such algorithms can be computationally intensive, particularly when sculpting 3D geometries . A field that is known to closely approximate the desired solution while still satisfying Maxwell’s equations could serve as a better target for the GS algorithm.
GS algorithms are also limited by the fact that they only shape the light intensity and not the phase or polarization of the incident field . Control of both angular and linear momentum is desirable for many optical tweezer applications and control of phase and polarization distributions could be implemented here through appropriate modification of the objective function. Although our method does not address how to physically create the desired fields, it could help to determine whether such a field is possible and whether it should be pursued by more intensive optimization of a DOE using, for example, a direct search algorithm .
Here we have considered particles that are well within the Rayleigh regime and are not expected to exert a significant scattering influence on the impinging field. For larger particles, this assumption will no longer hold and the desired field could be significantly altered if it has been designed assuming minimal interaction. Therefore, it may be advantageous to design fields that account for the presence of significant scatter. This could be accomplished within the present framework by making use of a T-matrix formulation and the translation and rotation properties of the VSWFs [43–46]. Using this procedure, the force on particles with a ≈ λ could be determined throughout the optimization volume. In addition, the T-matrix framework can be naturally extended to non-spherical particles.
Related work that has recently come to our attention combines the Debye-Wolf theory of light propagation and the Lorenz-Mie theory of light scattering to calculate the forces and torques present on particles in a holographic optical trap . This theory is fully vectorial and allows intensity, phase, and polarization effects to be modeled as a function of the phase and amplitude of the pixel elements governing the hologram. The authors report on the forces and torques in an optical vortex, a holographic ring trap, and a holographic line trap and comment on the untapped potential of “polarization engineering” for “control of highly structured torque fields”. While optimization routines and iterative algorithms have been used to generate tailored DOEs for specified light distributions, we are unaware of any work that seeks to generate tailored force fields. As we have shown here, such fields (including both linear and angular momentum) could be specified a priori and an evolutionary algorithm could be used to search the space of pixel amplitudes and phases governing the model given in  for those settings that most closely generate the desired field. Searching the space of VSWF coefficients using the procedure presented here might be a less expensive initial step prior to a full evolutionary algorithm optimization of the DOE.
In addition to tailored force fields for improved traps, we envision the possibility of creating “light wires” that could be used to guide particles through prescribed paths in a 3D volume of interest. Particles that start within specific regions would travel along designated light paths determined by the incident field. Such paths might be more complicated than the straight path generated by a donut beam. For example, the existence of spiral-type beams  bodes well for the existence of helical light wires. Such wires would differ from current 3D particle manipulation with optical tweezers because no active control of the incident beam would be required. This capability could be useful for optical construction of microdevices or fluidic sorting.
We have shown that Maxwell’s equations support the existence of a light field that approaches a desired light geometry: an optical tunnel. A vector spherical wavefunction basis was used to represent the incident field and an evolutionary algorithm (differential evolution) was used to synthesize the desired field in a least-squares sense. While we have not demonstrated how to physically produce the field, the presented framework might be used to find a target field for which a DOE could be created. Operating directly on the VSWF coefficients also allows extension of the procedure to particles whose size is on the order of a wavelength. Future work will seek to generate force fields capable of guiding particles through predefined 3D paths with no active control of the field.
References and links
3. A. Jonáš and P. Zemánek, “Light at work: the use of optical forces for particle manipulation, sorting, and analysis,” Electrophoresis 29, 4813–4851 (2008). [CrossRef]
5. D. L. Andrews, Structured Light and Its Applications: An Introduction to Phase-Structured Beams and Nanoscale Optical Forces (Elsevier, 2008). [PubMed]
6. P. J. Rodrigo, V. R. Daria, and J. Glückstad, “Four-dimensional optical manipulation of colloidal particles,” App. Phys. Lett. 86, 074103 (2005). [CrossRef]
7. T. Čižmár, V. Garcés-Chávez, K. Dholakia, and P. Zemánek, “Optical conveyor belt for delivery of submicron objects,” App. Phys. Lett. 86, 174101 (2005). [CrossRef]
8. V. Garcés-Chávez, K. Dholakia, and G. C. Spalding, “Extended-area optically induced organization of microparticies on a surface,” App. Phys. Lett. 86, 031106 (2005). [CrossRef]
9. J. Curtis, B. A. Koss, and D. G. Grier, “Dynamic holographic optical tweezers,” Opt. Commun. 207, 169–175 (2002). [CrossRef]
10. J. B. Wills, J. R. Butler, J. Palmer, and J. P. Reid, “Using optical landscapes to control, direct, and isolate aerosol particles,” Phys. Chem. Chem. Phys. 11, 8015–8020 (2009). [CrossRef] [PubMed]
11. J. Liesner, M. Reicherter, T. Haist, and H. J. Tiziani, “Multi-functional optical tweezers using computer-generated holograms,” Opt. Commun. 185, 77–82 (2000). [CrossRef]
12. J. Leach, G. Sinclair, P. Jordan, J. Courtial, M. Padgett, J. Cooper, and Z. Laczik, “3D manipulation of particles into crystal structures using holographic optical tweezers,” Opt. Express 12, 220–226 (2004). [CrossRef] [PubMed]
13. T. Čižmár, L. C. Dávila Romero, K. Dholakia, and D. L. Andrews, “Multiple optical trapping and binding: new routes to self-assembly,” J. Phys. B 43, 102001 (2010). [CrossRef]
15. J. Arlt, V. Garcés-Chávez, W. Sibbett, and K. Dholakia, “Optical micromanipulation using a Bessel light beam,” Opt. Commun. 197, 239–245 (2001). [CrossRef]
16. V. Garcés-Chávez, D. McGloin, W. Sibbett, and K. Dholakia, “Simultaneous micromanipulation in multiple planes using a self-reconstructing light beam,” Nature 419, 145–147 (2002). [CrossRef] [PubMed]
17. J. Arlt and M. J. Padgett, “Generation of a beam with a dark focus surrounded by regions of higher intensity: the optical bottle beam,” Opt. Lett. 25, 191–193 (2000). [CrossRef]
18. A. E. Chiou, W. Wang, G. J. Sonek, J. Hong, and M. W. Berns, “Interferometric optical tweezers,” Opt. Commun. 133, 7–10 (1997). [CrossRef]
19. P. C. Morgensen and J. Glückstad, “Dynamic array generation and pattern formation for optical tweezers,” Opt. Commun. 175, 75–81 (2000). [CrossRef]
20. E. R. Dufresne and D. G. Grier, “Optical tweezer arrays and optical substrates created with diffractive optics,” Rev. Sci. Instrum. 69, 1974–1977 (1998). [CrossRef]
21. M. Reicherter, T. Haist, E. U. Wagemann, and H. J. Tiziani, “Optical particle trapping with computer-generated holograms written on a liquid-crystal display,” Opt. Lett. 24, 608–610 (1999). [CrossRef]
22. K. Sasaki, M. Koshioka, H. Misawa, N. Kitamura, and H. Masuhara, “Optical trapping of a metal-particle and a water droplet by a scanning laser-beam,” App. Phys. Lett. 60, 807–809 (1992). [CrossRef]
23. C. Mio, T. Gong, A. Terray, and D. W. M. Marr, “Design of a scanning laser optical trap for multiparticle manipulation,” Rev. Sci. Instrum. 71, 2196–2200 (2000). [CrossRef]
24. K. Visscher, S. P. Gross, and S. M. Block, “Construction of multiple-beam optical traps with nanometer-resolution position sensing,” IEEE J. Sel. Top. Quantum Electron. 2, 1066–1075 (1996). [CrossRef]
25. M. T. Valentine, N. R. Guydosh, B. Gutiérrez-Medina, A. N. Fehr, J. O. Andreasson, and S. M. Block, “Precision steering of an optical trap by electro-optic deflection,” Opt. Commun. 33, 599–601 (2008).
26. R. Piestun, B. Spektor, and J. Shamir, “Wave fields in three dimensions: analysis and synthesis,” J. Opt. Soc. Am. A 13, 1837–1848 (1996). [CrossRef]
27. R. Piestun, B. Spektor, and J. Shamir, “Unconventional light distributions in three-dimensional domains,” J. Mod. Opt. 43, 1495–1507 (1996). [CrossRef]
28. G. Shabtay, “Three-dimensional beam forming and Ewald’s surfaces,” Opt. Commun. 226, 33–37 (2003). [CrossRef]
29. G. Whyte and J. Courtial, “Experimental demonstration of holographic three-dimensional light shaping using a Gerchberg–Saxton algorithm,” New J. Phys. 7, 117 (2005). [CrossRef]
30. O. Moine and B. Stout, “Optical force calculations in arbitrary beams by use of the vector addition theorem,” J. Opt. Soc. Am. B 22, 1620–1631 (2005). [CrossRef]
31. J. D. Jackson, Classical Electrodynamics (Wiley, 1999).
32. L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing (Wiley, 1985).
33. T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Heckenberg, “Multipole expansion of strongly focussed laser beams,” J. Quant. Spectrosc. Radiat. Transf. 79–80, 1005–1017 (2003). [CrossRef]
34. Y. Harada and T. Asakura, “Radiation forces on a dielectric sphere in the Rayleigh scattering regime,” Opt. Commun. 124, 529–541 (1996). [CrossRef]
35. S. H. Simpson and S. Hanna, “Rotation of absorbing spheres in Laguerre–Gaussian beams,” J. Opt. Soc. Am. A 26, 173–183 (2009). [CrossRef]
36. G. Gouesebet, J. A. Locke, and G. Gréhan, “Partial-wave representations of laser beams for use in light-scattering calculations,” Appl. Opt. 34, 2133–2143 (1995). [CrossRef]
37. L. W. Davis, “Theory of electromagnetic beams,” Phys. Rev. A 19, 1177–1179 (1979). [CrossRef]
38. R. Storn and R. Price, “Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces,” J. Global Optim. 11, 341–359 (1997). [CrossRef]
40. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).
41. G. C. Spalding, J. Courtial, and R. Di Leonardo “Holographic Optical Tweezers,” in Structured Light and Its Applications, D. L. Andrews, ed. (Elsevier, 2008), Chap. 6. [CrossRef]
43. C. H. Choi, J. Ivanic, M. S. Gordon, and K. Ruedenberg, “Rapid and stable determination of rotation matrices between spherical harmonics by direct recursion,” J. Chem. Phys. 111, 8825–8831 (1999). [CrossRef]
44. G. Videen, “Light Scattering from a Sphere Near a Plane Surface,” in Light Scattering from Microstructures, Lecture Notes in Physics Volume 534, F. Moreno and F. González, eds. (Springer, 2000), Chapter 5. [CrossRef]
45. T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Heckenberg, “Calculation of the T-matrix: general considerations and application of the point-matching method,” J. Quant. Spectrosc. Radiat. Transf. 79–80, 1019–1029 (2003). [CrossRef]
46. T. A. Nieminen, V. L. Y. Loke, A. B. Stilgoe, G. Knöner, A. M. Brańczyk, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Optical tweezers computational toolbox,” J. Opt. A 9, S196–S203 (2007). [CrossRef]