Tailored optical force fields using evolutionary algorithms

Colin C. Olson; Ross T. Schermer; Frank Bucholtz

doi:10.1364/OE.19.018543

1. Introduction

Optical trapping of particles using tightly focused laser beams was first demonstrated by Ashkin et al. [1] 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 [5]. 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 [13]. 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 [14], 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 [17]. The means by which these beams and traps are generated is varied, and includes multi-beam interference [18], the generalized phase contrast method [19], diffractive optical elements (DOEs) [20], holograms [21], 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 by

E (r, t) = R e [E (r) exp (- i ω t)]

where the complex amplitude E(r) satisfies

\nabla^{2} E (r) - k^{2} E (r) = 0

with k = n_mω/c, n_m the refractive index of the medium, ω the optical angular frequency, and c the speed of light in vacuum. We define a polar spherical coordinate system centered at the Cartesian coordinates x = 0, y = 0, z = 0 with position vector r = (x, y, z) = (r, θ, ϕ) and respective unit vectors x ̂ŷ,ẑ and r̂, $\hat{θ}$ , $\hat{ϕ}$ . The co-latitude, θ, is measured from the +z-axis and the azimuth, ϕ, is measured from the +x-axis toward the +y-axis. The VSWFs are a complete basis set of solutions to Eq. (2), thus, any arbitrary complex amplitude can be written as a linear combination of the VSWFs as follows:

E (r) = A \sum_{n = 1}^{\infty} \sum_{m = - n}^{n} [a_{n m}^{(1)} R g M_{n m} (k r, θ, φ) + a_{n m}^{(2)} R g N_{n m} (k r, θ, φ)]

where A is an overall field amplitude coefficient with dimensions of electric field and the expansion coefficients

a_{n m}^{(1)}

and

a_{n m}^{(2)}

are dimensionless. The magnetic field is given in turn by

H (r) = A {(\frac{ε_{r} ε_{0}}{μ_{r} μ_{0}})}^{1 / 2} \sum_{n = 1}^{\infty} \sum_{m = - n}^{n} [a_{n m}^{(1)} R g N_{n m} (k r, θ, φ) + a_{n m}^{(2)} R g M_{n m} (k r, θ, φ)]

where ε ₀ and μ ₀ are the permittivity and permeability of free space, respectively, and ε_r and μ_r are the dimensionless relative permittivity and permeability of the medium [30].

The singularity-free regular VSWFs are given by

R g M_{n m} (k r, θ, φ) = N_{n} j_{n} (k r) C_{n m} (θ, φ)

R g N_{n m} (k r, θ, φ) = \frac{j_{n} (k r)}{k r N_{n}} P_{n m} (θ, φ) + N_{n} (j_{n - 1} (k r) - \frac{n j_{n} (k r)}{k r}) B_{n m} (θ, φ)

where j_n(kr) is a spherical Bessel function and

N_{n} = 1 / \sqrt{n (n + 1)}

. The vector spherical harmonics are given by

B_{n m} (θ, φ) = \hat{θ} \frac{\partial}{\partial θ} Y_{n}^{m} (θ, φ) + \hat{φ} \frac{i m}{sin θ} Y_{n}^{m} (θ, φ)

C_{n m} (θ, φ) = \hat{θ} \frac{i m}{sin θ} Y_{n}^{m} (θ, φ) - \hat{φ} \frac{\partial}{\partial θ} Y_{n}^{m} (θ, φ)

P_{n m} (θ, φ) = \hat{r} Y_{n}^{m} (θ, φ)

where

Y_{n}^{m} (θ, φ)

is the normalized scalar spherical harmonic given by

Y_{n}^{m} (θ, φ) = \sqrt{\frac{2 n + 1}{4 π} \frac{(n - m)!}{(n + m)!}} P_{n}^{m} (cos θ) e^{i m φ}

and

P_{n}^{m} (\cdot)

is the associated Legendre function [30–33].

Choosing N_max 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 $a_{n m}^{(1)}$ and $a_{n m}^{(2)}$ (for notational convenience, we combine the $N_{m a x}^{2} + 2 N_{m a x}$ expansion coefficients into a single complex vector, $\tilde{a} = [a_{n m}^{(1)}, a_{n m}^{(2)}]$ ). Here we are not concerned with analysis, that is, determining the expansion coefficients given a defined field E(r) as discussed in [33]. 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 [34], the instantaneous gradient force is given by

F_{g r a d} (r, t) = 4 π n_{m}^{2} ε_{0} R^{3} (\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}}) \frac{1}{2} \nabla E^{2} (r, t)

where n_p is the index of the particle and R is its radius. Using the relation,

〈 E^{2} (r, t) 〉 = \frac{1}{2} {| E (r) |}^{2}

, where the angle brackets represent a time-average, the steady-state force is given by

\begin{array}{l} F_{g r a d} (r) & = & 〈 F_{g r a d} (r, t) 〉 \\ = & π n_{m}^{2} ε_{0} R^{3} (\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}}) \nabla {| E (r) |}^{2} \end{array}

where the gradient field, ∇|E(r)|², can be calculated numerically using finite differences if E(r) is known. The Poynting vector represents the instantaneous energy flux crossing a unit area per unit time and is given as

S (r, t) \equiv E (r, t) \times H (r, t)

with the time-dependent magnetic field given by

H (r, t) = R e [H (r) exp (- i ω t)] .

The scattering force on a particle is given by

F_{s c a t} (r) = \frac{n_{m}}{c} C_{s c a t} 〈 S (r, t) 〉

where

〈 S (r, t) 〉 = \frac{1}{2} R e [E^{*} (r) \times H (r)]

is the time-average of the Poynting vector and

C_{s c a t} = \frac{8}{3} π k^{4} R^{6} {(\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}})}^{2} .

is the scattering cross-section of the particle in the Rayleigh regime. Finally, the total force on the particle is simply

F_{t o t a l} (r) = F_{g r a d} (r) + F_{s c a t} (r) .

Thus, the force that a Rayleigh particle will encounter anywhere in a beam defined by a given set of expansion coefficients, ã, can be determined by numerically calculating the gradient field and using Eqs. (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, P_i, should be calculated for beams expanded as VSWFs. Moine and Stout [30] were careful to ensure proper normalization of beam power for particles with radii on the order of beam wavelength and questioned the idea that P_i ∝ ||ã||², where the double-bars represent the Euclidean norm. This proportionality was later shown by Simpson and Hanna [35] to be true for a wide variety of laboratory beams. By properly accounting for A in Eq. (21) of [35] we calculate the prefactor as

A^{2} = \frac{8 P_{i} k^{2}}{| | \tilde{a} | |} \sqrt{\frac{μ_{m}}{ε_{m}}} .

where μ_m = μ_rμ ₀ and ε_m = ε_rε ₀ are the permeability and permittivity of the medium.

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 [34]. The VSWF coefficients required to represent the Gaussian beam are found by using the development presented in [30]. First, the analytic expressions for the plane-wave coefficients are given as

\begin{array}{l} p_{n, 1}^{(1)} & = & i^{n} \sqrt{π (2 n + 1)} (i \hat{θ} + \hat{φ}) \cdot {\hat{e}}_{i} \\ p_{n, 1}^{(2)} & = & i^{n} \sqrt{π (2 n + 1)} (i \hat{θ} + \hat{φ}) \cdot {\hat{e}}_{i} \\ p_{n, - 1}^{(1)} & = & i^{n} \sqrt{π (2 n + 1)} (i \hat{θ} - \hat{φ}) \cdot {\hat{e}}_{i} \\ p_{n, - 1}^{(2)} & = & i^{n} \sqrt{π (2 n + 1)} (- i \hat{θ} + \hat{φ}) \cdot {\hat{e}}_{i} \end{array}

where all coefficients with m ≠ |1| are zero, ê_i is a polarization vector with ||ê_i|| = 1, and ã has been replaced by

\tilde{p} = [p_{n m}^{(1)}, p_{n m}^{(2)}]

to explicitly illustrate that the coefficients represent a plane wave. The VSWF coefficients describing an axisymmetric tightly focused beam can be determined by

{\tilde{a}}_{n} = g_{n} {\tilde{p}}_{n}

where the g_n are beam-shape coefficients [36]. A Davis-type model of a Gaussian beam [37] allows for nonparaxial corrections by taking s ² power series approximations of the vector potential, where the dimensionless beam shape parameter is defined by

s \equiv \frac{1}{k w_{0}}

and w ₀ is the beam waist radius in the paraxial approximation. The Davis model is useful here because it allows the g_n to be calculated analytically by either a first-, third-, or fifth-order expansion given, respectively, as

\begin{array}{l} \begin{array}{l} g_{n}^{(1)} & = & exp [- s^{2} (n - 1) (n + 2)] \\ g_{n}^{(3)} & = & g_{n}^{(1)} + exp [- s^{2} (n - 1) (n + 2)] (n - 1) (n + 2) s^{4} [3 - (n - 1) (n + 2) s^{2}] \\ g_{n}^{(5)} & = & g_{n}^{(3)} + exp [- s^{2} (n - 1) (n + 2)] \end{array} \\ \times {(n - 1)}^{2} {(n + 2)}^{2} s^{8} [10 - 5 (n - 1) (n + 2) s^{2} + 0.5 {(n - 1)}^{2} {(n + 2)}^{2} s^{4}] . \end{array}

Given (20), (21), and (23) we have the requisite coefficients ã for a focused Gaussian beam and are now in a position to compare F _grad (r) and F _scat (r) calculated using (12) and (15) to analytical representations of the same.

The analytical description of the scattering force as calculated in [34] is given by

F_{s c a t, z}^{(a)} (r) = \hat{z} \frac{n_{m}}{c} \frac{8}{3} π k^{4} R^{6} {(\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}})}^{2} (\frac{2 P_{i}}{π w_{0}^{2}}) \frac{1}{1 + {(2 \tilde{z})}^{2}} exp [- \frac{2 ({\tilde{x}}^{2} + {\tilde{y}}^{2})}{1 + {(2 \tilde{z})}^{2}}] .

where

(\tilde{x}, \tilde{y}, \tilde{x}) = (x / w_{0}, y / w_{0}, z / k w_{0}^{2})

. The scattering force is only assumed to be in the ẑ-direction because the model breaks down with highly focused beams and the assumption of only one transverse component of the electric field is accurate given the value of s = 0.012 used in their work. The analytical descriptions of the gradient forces in the x̂-, ŷ-, and ẑ-directions are given, respectively, by

F_{g r a d, x}^{(a)} (r) = - \hat{x} \frac{2 π n_{m} R^{3}}{c} (\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}}) \frac{4 \tilde{x} / w_{0}}{1 + {(2 \tilde{z})}^{2}} (\frac{2 P_{i}}{π w_{0}^{2}}) \frac{1}{1 + {(2 \tilde{z})}^{2}} exp [- \frac{2 ({\tilde{x}}^{2} + {\tilde{y}}^{2})}{1 + {(2 \tilde{z})}^{2}}]

F_{g r a d, y}^{(a)} (r) = - \hat{y} \frac{2 π n_{m} R^{3}}{c} (\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}}) \frac{4 \tilde{y} / w_{0}}{1 + {(2 \tilde{z})}^{2}} (\frac{2 P_{i}}{π w_{0}^{2}}) \frac{1}{1 + {(2 \tilde{z})}^{2}} exp [- \frac{2 ({\tilde{x}}^{2} + {\tilde{y}}^{2})}{1 + {(2 \tilde{z})}^{2}}]

\begin{array}{l} F_{g r a d, z}^{(a)} (r) = - \hat{z} \frac{2 π n_{m} R^{3}}{c} (\frac{n_{p}^{2} - n_{m}^{2}}{n_{p}^{2} + 2 n_{m}^{2}}) \frac{8 \tilde{z} / (k w_{0}^{2})}{1 + {(2 \tilde{z})}^{2}} [1 - \frac{2 ({\tilde{x}}^{2} + {\tilde{y}}^{2})}{1 + {(2 \tilde{z})}^{2}}] \\ \times (\frac{2 P_{i}}{π w_{0}^{2}}) \frac{1}{1 + {(2 \tilde{z})}^{2}} exp [- \frac{2 ({\tilde{x}}^{2} + {\tilde{y}}^{2}}{1 + {(2 \tilde{z})}^{2}}] . \end{array}

where we have added the factor of 2 to (25) and (26) that was missing from Eqs. (17) and (18) of [34]. A moderately-focused, x-polarized (ê_i = [1 0]) beam acting on polystyrene spheres in water was assumed in [34] with physical parameters n_p = 1.592, n_m = 1.332, w ₀ = 5 μm, P_i = 100 mW, and a vacuum wavelength λ ₀ = 0.5145 μm. We use the same values and calculate the Gaussian fields using N_max = 130. Analytical and computational results for the Gaussian beam are compared in Figure 1.

Fig. 1 Comparison between computational results (N_max = 130) and analytical results derived in [34]. All fields are calculated at the z = 0 plane. Analytical gradient force in the (A) x-direction and (B) y-direction. (C) Analytical scattering force in the z-direction. Computationally determined gradient field in the (D) x-direction and (E) y-direction. (F) Computationally determined scattering force in the z-direction. Gradient and scattering forces are negligible in the directions that are not shown. Small errors that appear in the outermost regions of the computational results could be eliminated by increasing N_max.

Download Full Size | PDF

We begin building the target force field, F _T (r), by generating a Gaussian beam for a particle with R = 50 nm. Fewer coefficients (N_max = 12) are required to accurately represent a more focused beam with w ₀ = λ_m where λ_m = λ ₀/n_m. 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 ₌₀. 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.

Fig. 2 (A) Target force distribution, F _T (r), for the z = 0 plane. (B) Target force distribution repeated at several z planes. The displayed force vectors have been scaled for visualization.

Download Full Size | PDF

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 letting

α^{(1)} = [a_{1, - 1}^{(1)}, a_{1, 0}^{(1)}, a_{1, 1}^{(1)}, a_{2, - 2}^{(1)}, a_{2, - 1}^{(1)}, \dots, a_{N_{m a x}, N_{m a x - 1}}^{(1)}, a_{N_{m a x}, N_{m a x}}^{(1)}]

α^{(2)} = [a_{1, - 1}^{(2)}, a_{1, 0}^{(2)}, a_{1, 1}^{(2)}, a_{2, - 2}^{(2)}, a_{2, - 1}^{(2)}, \dots, a_{N_{m a x}, N_{m a x - 1}}^{(2)}, a_{N_{m a x}, N_{m a x}}^{(2)}]

and defining

v = [R e {α^{(1)}}, I m {α^{(1)}}, R e {α^{(2)}}, I m {α^{(2)}}]

If we seek to match all force components comprising F _T (r) over the sampled volume of interest, then we are solving an optimization problem given by

min_{v} \sum_{r_{n} \in R} {(F_{t} (r_{n}; v_{t}) - F_{T} (r_{n}))}^{2}

where R is the set of all discretized points, r _n, that we have sampled in the volume of interest and F _t (r) is a test force field associated with a given set of expansion coefficients, v _t. Designating 𝒜 as the number of coefficients corresponding to a given N_max, the dimension of the search space for the optimization (31) is 2𝒜 which can become quite large in the fully general case of

𝒜 = N_{m a x}^{2} + 2 N_{m a x}

. If the space is restricted to axisymmetric solutions (which we impose for this study), then only the coefficients with |m| = 1 are required and the number of coefficients scales as 𝒜 = 4N_max [33].

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 N_max 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) [38]. 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 F_s. 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 F_s = 0.9 and CR = 0.5 as recommended in [38], 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.

Fig. 3 (A) The best discovered force distribution shown for the z = 0 plane. (B) A zoom of (A). The displayed force vectors have been scaled for visualization.

Download Full Size | PDF

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.

Fig. 4 Particle trajectories where all trajectories have been initiated from the z/λ = −1 plane. The first column is a three-quarter view of the trajectories with three slices of the force field included. The second column is a projection along the positive z-axis with the force field at z/λ = −1 superimposed. All figures show particle trajectories beginning from x/λ = −2 (x = −1.029 μm), x/λ = −1 (x = −0.5145 μm), and x = 0. The first row (AB) corresponds to all initial coordinates with y/λ = −1.5 (y = −0.7717 μm). The second row (C-D) corresponds to initial coordinates with y = 0. The initial point is represented by a circle and the final point by a square. The triangle is a central point and indicates the direction of particle movement.

Download Full Size | PDF

Fig. 5 Particle trajectories where all trajectories have been initiated from the z/λ = −1 plane. Each column has the same initial x-coordinates for all trajectories. Trajectories in the first column start at x/λ = −2 (x = −1.029 μm), the second at x/λ = −1 (x = −0.5145 μm), and the third at x = 0. The first row (A-C) corresponds to all initial coordinates with y/λ = −1.5 (y = −0.7717 μm). The second row (D-F) corresponds to initial coordinates with y/λ = −1.0 (y = −0.5145 μm). The third and fourth rows have initial conditions with y/λ = −0.5 (y = −0.2572 μm) and y = 0, respectively. Three slices of the force field are included. The initial point is represented by a red circle and the final point by a red square. The red triangle is a central point and indicates the direction of particle movement.

Download Full Size | PDF

Fig. 6 (A)–(F) show the z-component of the Poynting vector at the z/λ = −2, z/λ = −1.4, z/λ = −0.8, z/λ = 0.8, z/λ = 1.4, and z/λ = 2 planes, respectively. (G)–(I) show, respectively, the x-, y-, and z-components of the Poynting vector at the z = 0 plane. (A) and (F) are the only plots of S_z with negative values.

Download Full Size | PDF

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 [14] and Bessel beams [15]. 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 [39], however, it is generally not known a priori whether the desired field is consistent with Maxwell’s equations [26]. The algorithms used to produce the DOE, such as Gerchberg-Saxton (GS) [40], 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 [41]. 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 [29]. 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 [42].

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 [47]. 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 [47] 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 [27] 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.

6. Conclusion

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

1. 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]

2. D. J. Stevenson, F. Gunn-Moore, and K. Dholakia, “Light forces the pace: optical manipulation for biophotonics,” J. Biomed. Opt. 15, 041503 (2010). [CrossRef] [PubMed]

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]

4. D. G. Grier, “A revolution in optical manipulation,” Nature 424, 810–816 (2003). [CrossRef] [PubMed]

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]

14. K. T. Gahagan and G. A. Swartzlander Jr., “Optical vortex trapping of particles,” Opt. Lett. 21, 827–829 (1996). [CrossRef] [PubMed]

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]

39. R. Piestun and J. Shamir, “Control of wave-front propagation with diffractive elements,” Opt. Lett. 19, 771–773 (1994). [CrossRef] [PubMed]

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]

42. M. A. Seldowitz, J. P. Allebach, and D. W. Sweeney, “Synthesis of digital holograms by direct binary search,” Appl. Opt. 26, 2788–2798 (1987). [CrossRef] [PubMed]

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]

47. B. Sun, Y. Roichman, and D. G. Grier, “Theory of holographic optical trapping,” Opt. Express 16, 15765–15776 (2008). [CrossRef] [PubMed]

Tailored optical force fields using evolutionary algorithms

Abstract

1. Introduction

2. Field expansion

3. Force field calculation

4. Design procedure

5. Results and discussion

6. Conclusion

References and links

Cited By

Figures (6)

Equations (31)

Optics Express