1. Introduction

Trapping of particles has been of great interest to science and engineering [1] since Ashkin demonstrated trapping of dielectric particles using radiation forces produced by continuous laser beams [2]. A simple method to trap micrometer-sized dielectric particles is performed using dual counter-propagating divergent beams, which enables trapping of single particles [3,4], and self-organized multi-particle arrays [5-11]. The self-organized multi-particle trap offers an interesting system for studying optical interactions between particles, with applications towards contact-free storage of biological cells [6] and ion traps for quantum computing [5,12]. Recently, there has been much interest in understanding the underlying physics behind optical binding in multi-particle traps [5-11].

A demonstration of self-organization in an dual-beam optical trap was carried out by Tatarkova et al. using counter-propagating laser beams [5]. They presented experimental results of equilibrium positions for up to seven trapped particles in an array and discussed modes of oscillation, which are similar to those involving ion traps. Later, Singer et al. undertook experiments containing different particle sizes to investigate the dependency of equilibrium spacing on the particle size [6]. They determined that the particles were equidistantly separated and showed a strong dependency of equilibrium spacing on particle size. Also, observed was a slight dependency of particle spacing on the refractive index difference between the particle and surrounding medium, on fiber spacing, and on the number of trapped particles. Using a model based on first-order plane wave scattering from point scatterers inside an infinite one-dimensional array, they predicted equilibrium spacing as a function of particle size in a size regime around the used wavelength (0.5λ<D<2 λ, D is the trapped particle diameter). As this model used an infinite array of point scatterers, several important phenomena could not be estimated, such as the dependency of equilibrium position on the number of particles, dependency on particle-host medium refractive index, and dependency on opposing fiber separation. Later, McGloin et al. proposed a numerical model based on a two-component approach, involving gradient and scattering forces, combined with the paraxial wave equation, where they calculated equilibrium spacing for a two- and three-particle system [7]. They showed numerically the occurrence of inhomogeneous equilibrium spacing as a function of the number of particles, particularly that the spacing for a three-particle array is smaller than that for two-particle array. Since their model relies on the paraxial wave theory, it fails for small particles and for a high refractive index contrast between the particle and surrounding medium. Using a similar model, Metzger et al. simulated bistability and hysteresis in optical binding between two particles [8], and further analyzed restoring forces on two optically bound particles located around their equilibrium positions while also including hydrodynamic interactions (HDI) [11]. In addition, a coupled diode method has also been applied to a 2-sphere array where several important phenomena were demonstrated, namely, the sensitivity of equilibrium positions on refractive index contrast and the effects of standing waves caused by the incident beam and scattered fields [13].

Although several studies have contributed to understanding the behavior of trapped particles in a dual-beam system, further investigations are still necessary, to explain the phenomena associated with large number multi-particle arrays [14, 15]. Currently, no theory has explained fully the occurrence of inhomogeneous particle spacing, both for a particle number dependency and a dependence on inter-array particle positions (i.e. inner and outer inter-particle spacings differ for a fixed number particle array), and the spontaneous onset of oscillations observed in the dual beam trap [5].

In this work, trapping is demonstrated in an integrated microfluidic-fiber optic environment to study trapping physics for a large number-particle array. Two types of inhomogeneous properties of equilibrium spacing are observed; one is that, for a fixed number particle-array, the inner inter-particle spacing is smaller than the outer inter-particle spacing, the other is that the averaged spacing of an I-particle array decreases as I increases. Furthermore, the theory describes the onset of spontaneous oscillation for a critical number of particles, which has not been theoretically analyzed for a multi-particle array. To analyze the experimentally observed phenomena, optical forces acting on the trapped particles are numerically calculated using the Maxwell stress tensor with Generalized Multipole Technique (MST-GMT) [16, 17]. Since the fields inside and outside the particles in the GMT are expanded by spherical vector wave functions, light scattering in the Rayleigh and Mie regimes can be accurately calculated while including a high refractive index contrast. The MST-GMT results are verified using FDTD, which is less efficient for this system. Using the MST-GMT, dynamic simulations are performed which include hydrodynamic interactions between particles [18]. The influence of hydrodynamic interaction is demonstrated by setting these interactions to zero and comparing the results. The numerical results show inhomogeneous spacing between trapped particles that agrees well quantitatively with experimental results. Also, the dynamic simulation shows that the critical onset of oscillation is caused by a misalignment of the beam axes, as is found in the experiment.

2. Experiment

2.1 Chip fabrication and experimental setup

A Polydimethysiloxane (PDMS) optofluidic chip was fabricated with two integrated, counter-directional optical fibers. The chip was patterned using CAD software and fabricated using established soft-lithography techniques [19]. Following chip assembly, optical fibers were inserted into two opposing 250×250 µm² vias where they come in contact with the perpendicular microfluidic channel used for particle delivery. The axial separation of the two cleaved fibers was set at 160 µm, whereas the layout and molding of the fiber vias ensured radial pre-alignment. Each fiber was coupled to a separate laser diode source, both operating at a wavelength of 980 nm and emitting 100 mW.

The polystyrene spheres (Bangs Laboratory Inc.) of 1-micron-diameter were purchased in solution and so the refractive index does not change appreciably over the duration of the experiment. The refractive index of a polystyrene sphere in water suspension was previously found to be 1.574 [20]. The spheres were fluorescein-labeled to enhance imaging. The particle concentration in water (n=1.33) was approximately 0.01% and the particle-solution delivery was controlled with a standard 1 mL syringe press fitted to polymer tubing.

2.2 Measurement and imaging

Trapped stationary and oscillating particles were imaged with a fluorescence microscope containing a 10×(NA=0.30) objective (DMLB, Leica Microsystems, Germany). Still images at 30fps were acquired from a CCD camera coupled to the microscope objective. Inter-particle spacing measurements and oscillation amplitudes were determined using known size calibrations and pixel counting. Similarly, oscillation periods were measured from known image capture rates.

2.3 Inhomogeneous spacing and oscillation

Figure 1 shows images of particle array trapping for 10 particles, and dynamic oscillation for 13 particles. A video is also provided which clearly shows the onset of large-amplitude oscillations after a critical number of particles is reached (in this case, 13 particles). The outerparticle oscillation amplitude was measured to be 20 µm and the period was 0.5 seconds.

 

Fig. 1. Fluorescence microscopy images of a stably trapped 10-particle array containing 1-micron fluorescent polystyrene spheres and a 13-particle oscillating array. A video demonstrating the transition from a 10- to 13-particle array is also provided.

Download Full Size | PDF

3. Optohydrodynamic theory

To analyze the experimentally obtained results, an optohydrodynamic numerical model is provided. Fig. 2 shows polystyrene (PS) spheres of radius a trapped between opposing fibers separated by a distance D _bf. The divergent beams (DB1 and DB2) from the optical fibers (OF1 and OF2) are modeled by a Gaussian beam. To make the model realistic, we consider the case of transversely misaligned beam axes (D _off≠0), in addition to the case of perfectly aligned beam axes (D _off=0), where D _off is the distance between the beam axes (BA1 and BA2). As will be shown later, this misalignment results in spontaneous oscillation of the particle array.

 

Fig. 2. Schematic used to model of the dual-beam fiber trap.

Download Full Size | PDF

3.1 Optical force calculation

To calculate the inhomogeneous spacing and the spontaneous oscillation observed in the experiment, the optical force acting on particles in the dual-beam trap needs to be calculated. The time averaged optical force F _i acting on the i ^th particle can be calculated by integrating the Maxwell stress tensor T⃡ over a closed surface S_{M, i} which surrounds the particle [21]:

where v _M is a unit normal vector of the closed surface. The Maxwell stress tensor (MST), T⃡ is given by

where ε and µ are the permittivity and the permeability of the surrounding medium, respectively, * represents the complex conjugate of, EE* is a dyadic product, I⃡ is the unit tensor, and Re denotes the real part of. The electric and magnetic fields E and H are, respectively, given by

where E ⁱⁿ and H ⁱⁿ are the incident fields, i.e. divergent beams from the fiber facets, and E ^sc and H ^sc are the scattered field by particles.

In order to obtain the scattered fields E ^sc and H ^sc numerically, we apply the Generalized Multipole Technique (GMT) [16, 17] to the multi-sphere system. To apply this method, we first express the approximate scattered electric field at a point r, E ^sc _N(r), as

where I is the total number of particles, r _i denotes a position r in the ith local coordinate system with origin O_i located at the centre of the i ^th particle, and E ^sc _i,N(r _i), is the scattered electric field from the i ^th particle. Following the GMT, we expand E ^sc _i,N(r _i) as

where m ⁽⁴⁾ _nm and n ⁽⁴⁾ _nm are the Hankel function-based spherical vector wave functions (SVWFs) [22], N is a truncation size of the expansion that is given depending on a desired accuracy of the solution, and aⁱ_nm and bⁱ_nm are unknown coefficients to be determined. To calculate these coefficients, we also need to expand the transmitted field inside the i ^th particle, E ^tr _i,N(r _i), by SVWFs as

where m ⁽¹⁾ _nm and n ⁽¹⁾ _nm are Bessel function based SVWFs [22], and cⁱ_nm and dⁱ_nm are unknown coefficients. Using Faraday’s law and the relation, ∇×m=k n and ∇×n=k m, where k is the wave number, the scattered magnetic field from the i ^th particle, H ^sc _i,N(r _i), and the transmitted magnetic field inside the i ^th particle, H ^tr _i,N(r _i), are given by

where Z_out and Z_in are the intrinsic impedance of surrounding medium and that inside the i ^th particle, respectively.

The unknown coefficients, aⁱ_nm, bⁱ_nm, cⁱ_nm, and dⁱ_nm are determined by matching the boundary conditions on the surface of the particles. Here we numerically match the boundary conditions in a least squares sense [17], i.e. we calculate the unknown coefficients that minimize a squared norm as

where v _P is unit normal vector on the particle surface S_{P, j}. Minimization of the squared norm results in a K×K (K=4IN(N+2)) matrix equation [17]. The discretization of the squared norm is important from a numerical point of view and the discretization rule can be found elsewhere [17]. Once we obtain the expansion coefficients by solving the matrix equation, we can determine the field distributions inside and outside the particles by substituting these coefficients into Eq. (5) - (8). For a multi-particle array, the convergence on the number of basis functions used in the expansion was checked, and the typical number required for our calculations was N=7-10. By integrating the MST on a closed surface S_M,i outside the particles numerically, we can obtain the optical force acting on the i ^th particle. As a closed surface S_M,i, we choose a concentric sphere surface with radius b (>a) which surrounds only the i ^th particle as shown in Fig. 2. In our simulation, we assume that the counter-propagating beams are mutually incoherent [6, 9]; this is a reasonable assumption given that the laser diodes are separated by over 1 m of fiber in the experiments.

3.2 Dynamic simulation

To calculate the equilibrium positions for a large-number particle array and to investigate the dynamics of the particles in the array, we perform a dynamic simulation. The dynamic simulation can be performed by solving the equation of motion,

where R _i and v _i are the center position and the velocity of the i ^th particle, respectively. In our simulation, we approximate the velocity v _i using the Oseen tensor H⃡ [18],

where F _j is the optical force acting on the j ^th particle and

Here ζ=6πηa is the friction coefficient, η is the viscosity of the surrounding medium, δ_ij is Kronecker delta, and r_ij=|r _ij|=|r _i-r _j|. To evaluate Eq. (10) numerically, we use the Euler method, i.e. the position at time t=(m+1)Δt is determined by

where Δt is an appropriate time increment.

4. Numerical results and discussion

In the following numerical calculations, we use a fifth-order Gaussian beam [23] to simulate the divergent beam emitted from the optical fibers. The waist radius and the wavelength are set at ω₀=3.0 µm and λ=980 µm in free space, respectively, and the distance between the two fibers is set at D _bf=160 µm. The radius and refractive index of the particles are 0.5 µm and 1.574, respectively, and the refractive index of the surrounding water is 1.33.

4.1 Equilibrium spacing of two- and three- particle arrays

The equilibrium spacing for the two- and three-particle array can be obtained by calculating the optical force as a function of distance between neighboring particles and locating a point where the optical forces are zero. To verify the proposed MST-GMT numerical method, equilibrium spacing for the two- and three-particle array are compared with experimental results and with numerical results obtained by MST using a finite difference time domain method (MST-FDTD). For cases requiring MST-FDTD, a commercial FDTD software package, FDTD Solutions, provided by Lumerical Solutions Inc. [24] was used. In this section, we consider the case of perfectly aligned beam axes (D _off≠0).

Fig. 3(a) shows optical forces acting on a pair of particles when the spacing between particles, d, is changed from 2 µm to 20 µm. A modulation was found in the force calculations, as reported previously [13]. For the conditions of our experiment, the modulation is small so it does not modify significantly the large-scale dynamics. As shown in Fig. 3(a), good agreement exists between the MST-GMT and MST-FDTD method. As the spacing, d, is increased, the absolute value of the forces acting on particles 1 and 2 decrease and the forces take opposite sign after crossing the zero force line around d=9.4 µm. For d<9.4, particles repel each other, while for d>9.4 µm, particles attract. For d=9.4 µm, the net force acting on the particles is zero, and the calculated equilibrium spacing for the two particles is d _2,cal=9.4 µm. This spacing is in good agreement with the experimental result, d _2,exp=9.6 µm. In Fig. 3(b), forces acting on a three-particle array are shown, where the centre particle is fixed at z=0 and the distance, d, between this particle and the neighboring particles is varied. In this case, the calculated equilibrium spacing is around d _3,cal=8.3 µm and is also in good agreement with experimental result d _3,exp=8.0 µm. In agreement with Ref. [5], our experimental and numerical results show inhomogeneous equilibrium spacing of a two- and a three-particle array, i.e., the equilibrium spacing of a two-particle array is larger than that of a three-particle array.

Equilibrium forces were calculated for different refractive index values (n=1.64, 1.74, 1.84), however, substantial changes in the results were not observed.

 

Fig. 3. Optical forces acting on trapped particles as a function of distance between neighboring particles. Solid lines and blue dots are results obtained using MST-GMT and MST-FDTD, respectively. (a) Two-particle array. (b) Three-particle array.

Download Full Size | PDF

4.2 Equilibrium positions of an N-particle array

To calculate the equilibrium positions of an array consisting of more than three-particles, the dynamic simulation proposed in Section 3.2 is used. The case of perfectly aligned beam axes (D _off=0) is considered, where particle motion is confined to the z-axis.

Fig. 4(a) shows the results of this dynamic simulation. Simulation results with and without hydrodynamic interaction are displayed simultaneously, where red lines represent the inclusion of HDI, and blue lines without. The force without HDI is calculated by neglecting the second term in Eq. (12). This does not correspond to a real physical system, but it is useful for observing quantitative changes in the dynamics that result from modifications to the physical model. At t=0, the 1^st particle is initially located at an equilibrium position, z _1,1=0, and a 2^nd particle is introduced 20 µm away from the 1^st particle, at z=20 µm, to model an incoming particle. As the 2^nd particle is introduced, the two particles shift in the negative zdirection and settle at their new equilibrium positions, z _2,1=-4.7 µm and z _2,2=4.7 µm, which yields an equilibrium spacing d _2,dyn=9.4 µm. This result is identical to results obtained in Section 4.1, where d _2,cal=9.4 µm. Once the two particles have settled, a 3^rd particle is introduced 20 µm from the 2^nd particle (z=4.7+20 µm) at t/η(T)/10³≈0.9 m²/N, where η(T) is the viscosity of water as a function of temperature T [25]. The two previously stable particles shift in the negative z-direction and settle to their new equilibrium positions, z _3,1=-8.3 µm, z _3,2=0 and z _3,3=8.3 µm, for the 1^st, 2^nd, and 3^rd particles, respectively, resulting in an equilibrium spacing d _3,dyn=8.3 µm. This result also matches the value obtained in Section 4.1, where d _3,cal=8.3 µm. Extending this procedure to a 13-particle array, which was the critical number in the experiments that lead to array oscillations, equilibrium positions of a 13-particle array were calculated as shown in Fig. 4(b). At t=0, 12 particles are initially located at their equilibrium positions and a 13^th particle is introduced 20 µm from the 12^th particle. As in the previous case, the 13 particles settle to their new equilibrium positions and the onset of the spontaneous oscillation does not occur for this model containing perfectly aligned beam axes. The results with and without HDI do not show significant deviation for a collinear dual-fiber arrangement; however, the results will be shown to differ significantly for offset beams. Even for perfect alignment, it is clear that the HDI has a stronger influence for more closely spaced particles, as expected from Eq. 12. For example, the deviation between the equilibrium-approach with and without HDI for 4 particles is greater than for just two as in Fig. 4(a).

Figure 5(a) shows that the averaged inter-particle spacing, d _av, decreases in a nonlinear way as the total number of trapped particles, I, increases. The calculated theoretical results agree well with those measured in the experiment. Fig. 5(b) shows that, for a given number of trapped particles in the array, the outer particle spacing, d _out, is larger than the inner spacing, d _inn, near the array centre, and that the difference between the outer and inner spacing increases with increasing particle numbers. Again, numerical and experimental results are in close agreement.

 

Fig. 4. Trajectories of particles obtained by the dynamic simulation. Red and blue lines are the results with and without HDI, respectively. (a) From two- to four- particle array. (b) Thirteenparticle array.

Download Full Size | PDF

 

Fig. 5. Inhomogeneous properties of inter-particle spacing. (a) Averaged inter-particle spacing d _av as a function of the number of particles I. A standard deviation plot is applied to the experimental results. (b) The inner most inter-particle spacing d _inn and the outer most inter-particle spacing d _out, as a function of I. Green and red positions denote theoretical and experimental results, respectively.

Download Full Size | PDF

4.3 Dynamic simulation of the spontaneous oscillation

In this section, we consider a transversely misaligned model (D _off≠0) and show that the spontaneous oscillation observed in the experiment can be caused by this beam misalignment. We assume that the axis of a beam propagating in the positive z-direction lies on the z-axis and that propagating in the negative z-direction is located on a parallel axis with an offset in the positive x-direction, D _off, as shown in Fig. 2. Due to this beam offset, the particles can move in the x-z plane.

Calculations in this section are performed using the procedure mentioned in the previous section. First, a particle is positioned at the equilibrium position, i.e. (x₁, z₁)=(D _off/2, 0), and a second particle is introduced at (x₂, z₂)=(D _off, 20 µm) to model an incoming particle. Once the second particle is added, the two particles move towards their new equilibrium positions and settle as described previously. After settling, a third particle is added and this procedure is iterated for up to 7 particles.

 

Fig. 6. Simulation results for a 7-particle array. (a) Introduction of a 7^th particle into a system where the beam axes are offset by 4 µm. The arrow indicates the trajectory of the 7^th particle. (b) Final equilibrium positions for the 7 particles after undergoing damped oscillation induced by the introduction of the 7^th particle. (c) The resulting self-sustained oscillation after increasing the beam offset to 5 µm. A video demonstrating these results is also provided. [Video 2][Video 3]

Download Full Size | PDF

The dynamic simulation for different offsets and particle numbers were performed. The general properties of the trajectories are as follows. For a small number of particles, the particle trajectories follow an over damped behavior, as in the case of the perfectly aligned model shown in Fig. 4. As the particle numbers increase, damping of the particle oscillation decreases, and for a critical number of particles, stable oscillations occur. Fig. 6 and the accompanying video show 7 particles in motion for transverse offsets of 4 µm and 5 µm. In case of a 4 µm offset, the oscillations are critically damped and a stable final arrangement results. For the 5 µm offset, self-sustained oscillations are clearly visible.

The origin of the self-sustained oscillation in a viscously over-damped media can be explained by considering the different components of optical force acting on each particle. In Fig. 6, the z-displacement corresponds to varying the potential energy of the particle from the optical scattering force [26]. Similarly, the x-displacement corresponds to varying the potential of a particle from the optical gradient potential, which increases with increasing distance from the beam center [26]. This allows for an interpretation of the observed oscillations: the potential energy is exchanged between the optical scattering and gradient potentials to move the particle in an orbit in the x-z plane. This is fundamentally different from other coupled-particle systems where inertia plays an important role; in those systems there is an exchange between potential and kinetic energies. Here, the usual kinetic energy contribution is replaced by the gradient potential of the optical trap. It is clear from the optohydrodynamic theory that transverse symmetry breaking is required to produce the observed self-sustaining oscillations. The symmetry breaking allows for the coupling exchange between the scattering potential and the gradient potential, as mediated through optical interaction.

It is important to distinguish this oscillation from oscillations that may appear in systems due to inertial mass; for example, as described in a recent theoretical contribution [27]. In the present system, inertia is effectively damped out immediately and so there is no kinetic-energy contribution to the dynamics.

Since the behaviors of the oscillations depend on the number of particles and on the offset, trajectories were studied for various offsets, D _off, as a function of the number I, and quantitative estimates were made by fitting the trajectories to the damping relation Acos(Ωt)exp(-Λt)+B. Results of the angular frequency Ω and the damping coefficient Λ as a function of I are shown in Fig. 7, together with the results without HDI for reference. The threshold for self-sustained oscillation occurs when Λ becomes negative, which can be seen for the 7 particle case with a 5 µm offset. It is also clear that the oscillation frequency and the damping coefficient are larger when hydrodynamic interaction is included.

In the above examples, the offset was set greater than the experimental tolerances to observe oscillations for a small number of particles, since the calculation time for many particles is prohibitive for the current single-processor implementation of the numerical model. Although systematic studies of greater than 7 particles have not yet been performed, the base mechanism presented here concludes that the occurrence of spontaneous oscillation is the result of asymmetry, for example, from misalignment of the beam axes.

Since the theory uses GMT to calculate the scattered field, this method can be applied in both the Rayleigh and Mie regimes, including a high refractive index contrast. As the particle size and/or the refractive index contrast become larger, higher orders of basis expansion are needed, which means a larger sized matrix equation should be solved. In our calculation, after verifying the convergence of the solutions, the truncation size N was set to 8. As an example of computation time, when scattering by 7-particles was calculated (Fig. 6, supplementary video), the memory requirements for the matrix calculation was 480 MB and the computation lasted 16 hours using a 1.6 GHz PowerPC 970 single processor.

 

Fig. 7. Oscillation frequency and damping as a function of a number of the particles I for various offsets D _off. Squares and crosses represent results for D _off=5 µm and 4 µm, respectively, and red and blue lines represent results with and without HDI, respectively. (a) Frequency Ω, (b) Damping coefficient Λ.

Download Full Size | PDF

5. Conclusion

A theory was presented to explain the dynamics of optical trapping of polystyrene spheres in a dual-fiber system and was found to agree well with experimental observations. Two types of inhomogeneous properties in equilibrium spacing were shown for a trapped particle array; one is that, for a fixed I-particle array, the inner inter-particle spacing is smaller than the outer spacing, and the other is a nonlinear decrease of the averaged equilibrium spacing as the number of trapped particles increases. In addition, distinct spontaneous oscillations occurred for a large number of particles. The theory reproduced the features of the experiments by using a Maxwell stress tensor with the GMT while including hydrodynamic interactions. It was shown that a transverse misalignment of beam axes is the mechanism responsible for spontaneous self-sustained oscillation in our dual-beam trap.