Optoelectronic Tweezers under Arbitrary Illumination Patterns: Theoretical Simulations and Comparison to Experiment

Photovoltaic tweezers are a promising tool to place and move particles on the surface of a photovoltaic material in a controlled way. To exploit this new technique it is necessary to accurately know the electric field created by a specific illumination on the surface of the crystal and above it. This paper describes a numerical algorithm to obtain this electric field generated by several relevant light patterns, and uses them to calculate the dielectrophoretic potential acting over neutral, polarizable particles in the proximity of the crystal. The results are compared to experiments carried out in LiNbOs with good overall agreement. References 1. A. Ashkin, "History of optical trapping and manipulation of small-neutral particle, atoms and molecules," IEEETrapping of dielectric particles with light-induced space-charge fields," Appl.Two-dimensional dielectrophoretic particle trapping in a hybrid crystal/pdms-system," Opt.graphic surface gratings in iron-doped lithium niobate," Appl.Optical trapping and manipulation of metallic micro/nanoparticles via photorefractive crystals," Opt.Patterning of dielectric nanoparticles using dielectrophoretic forces generated by ferroelectric polydomain films," J.Role of particle anisotropy and deposition method on the patterning of nano-objects by the photovoltaic effect in LiNb03," Opt.Optimization of particle trapping and patterning via photovoltaic tweezers: role of light modulation and particle size," J. Phys. D Appl.Steady-state photorefractive gratings in LiNb03 for strong light modulation depths," IEEE J.


Introduction
The fast development of new research fields such as nanotechonology or biomedicine has required a parallel improvement of tools for manipulation of micro-and nano-objects [1][2][3].
Among them, electrokinetical methods based in electrophoretic (EP) or dielectrophoretic (DEP) fields are very appropriate for parallel trapping and patterning of many particles.An emergent methodology of this type [4][5][6][7] uses light-induced space charge fields generated in photorefractive materials by inhomogeneous illumination [8].The space charge fields responsible for particle trapping are particularly large in PR crystals exhibiting a large photovoltaic effect such as Fe:LiNb03.The micro-or nano-particles are deposited on the crystal surface by electrophoretic or dielectrophoretic forces induced by the evanescent space charge fields.The procedure does not require electrodes and it is usually called photovoltaic or, more generally, optoelectronic tweezers.The number of experimental works showing the capabilities of optoelectronic tweezers has been rapidly increasing in the last few years [4][5][6][7][9][10][11][12].Different dielectric and metallic nanoparticles, neutral or charged, have been used to generate the particle patterns showing the universality of the method.Conversely, the theoretical description of the physics involved is at the beginning.One early paper [13] described particle trapping by DEP forces but appear in a different context (induced by electric fields appearing in ferroelectric domain structure).More recently, a theoretical approach already addressed to optoelectronic tweezers has been developed in two contributions [14,15], although it is essentially focused to periodic light illumination.In [14], Villarroel et al. present the theoretical framework, including the calculation of evanescent fields and dielectrophoretic forces for isotropic and anisotropic particles.The analysis is applied to single gaussian beam and sinusoidal illumination although, in order to obtain simple analytical expression, drastic approximations are applied and higher harmonics of the field are neglected, keeping only the fundamental and the first harmonic.
Presently we address the development of a completely general numerical tool able to obtain de electric field that acts over the particles surrounding the material.The interaction with this electric field can be calculated from the electronic structure of the particles and depends on its multipolar expansion.In this work we present the DEP potential acting on neutral, spherical and isotropic particles by the induction of a dipole moment.
The numerical tool allows to gain a complete picture of the physical phenomena showing the physics of the trapping process providing the two-dimensional profiles of the different magnitudes involved: space charge, evanescent electric field and DEP forces.

Theoretical model
Although the mathematical framework is, in essence, the same used in previous works [14,15], the analytical approach of these works forced the authors to make several approximations and the expressions obtained are valid only for a limited number of light configurations.In this section we will describe the calculation of i) the bulk and the evanescent space charge field and ii) the dielectrophoretic potential.

Bulk and evanescent space charge fields
The photovoltaic electric field and the current density generated in the material due to illumination can be calculated with the well known generaltabletions of the photorefractive effect [16].These equations assume one single type of free charge carrier and a single impurity center, having concentrations n and N respectively, which is a good approximation under moderate light intensities.So the rate equations governing the evolution of the system under arbitrary illumination I(r,t) are given by [8]: where Np> and NA are the donors and acceptors concentrations respectively, s is the photoionization cross section, 7 the recombination constant, ./thecurrent density (including the drift, diffusion and photovoltaic contributions), fi the mobility, E the total electric field (associated to that generated by space charge plus the external contribution, if any), D= iiksT/qthe diffusion coefficient, L pv the photovoltaic transport length and upv the unit vector in the direction of the polar axis.Note that, for any time, Np> and NA satisfy the relationship NA + NP> = N, which stands for the total center concentration and is independent of the presence of light.Besides the one center approximation, we have also ignored the thermal ionization of donors.This set of equations will be solved numerically as explain in section 3.

Electrical potential on the surface of the crystal
The electrical field generated inside the crystal leads to a force acting on the particles in the surroundings.In general, this force reaches its maximum value at the surface of the crystal and decreases exponentially with the distance.It acts upon particles, molecules, etc by means of different mechanisms depending on their electrical nature.In this work we will focus on the calculation of the attractive potential experienced by neutral particles or molecules immersed in a non-polar fluid.Small particles react to an electric field adapting their electronic cloud and moving accordingly.This interaction can be described in terms of the multipolar expansion of the electronic cloud of the particle: charge, dipole, quadrupole, and higher order multipoles.For electrically neutral particles, with no permanent multipoles, the strongest interaction with the space charge field of the crystal is due to the induced dipole on the particle.In this work we will assume dn dt SND dt J = electrically neutral, isotropic particles, which are small compared with the periodicity (if any) of the electric field, so the dipole approximation can be safely applied [17,18].
In this situation the force associated to the coupling energy is: where p is the induced dipolar moment [4,13], ep(0) + 2e M (0) being r the particle radius and a the scalar particle polarizability for an isotropic spherical particle.The values of the dielectric permitivities of the fluid (M) and the particle (F) correspond to static fields.As DEP force is equal to the gradient of a spatial function, it is a conservative force, and DEP potential, VDEP, can be defined:

Parameters used in the simulations
In order to obtain numerical values and to compare with experiments, the calculations have been done for an X-cut iron-doped lithium niobate crystal.Table 1 shows the values of the parameters used in the simulations [14].In this geometrical configuration of the crystal the Caxis is parallel to the surface.In our model the C-axis is parallel to the Y-axis and the Z-axis is perpendicular to it.We also assume the sample is thin enough to discard the absorption and consider the same illumination at any depth.Drift, diffusion and photovoltaic field contribute to the current density along the Y-axis, whereas only drift and diffusion contribute to current density along Z-axis.This leads to an homogeneous space charge field along the Z-axis.10 24 m-J photovoltaic transport length Lpy 5x 10-lu m

Numerical calculations
The calculation of the final DEP potential can be split in two independent steps: calculation of the steady state of the space charge distribution for a given illumination pattern I(x,y) and calculation of the DEP potential caused by this field on the surface of the crystal and in the outer nearby.

Steady state of the space charge Geld
The steady state of the space charge field is obtained by the numerical resolution of Eqs.(1), ( 2) and (3) using the values of Table 1 and the desired light intensity pattern.As the spatial charge field is homogeneous along the Z-axis, at this stage the crystal is treated as a surface (i.e.representing the crystal surface (see Fig. 1(a)).To simplify the calculation we make use of the so-called adiabatic approximation, that assumes that the relaxation time of charge carriers (electrons) is much shorter than the characteristic time scale for the variation of the electric field.E (i.e.recombination time, TR, much shorter than relaxation time, To).In this case dn/dt & 0 and consequently dn/dt = f(AN D , I, n{x,y)) = 0 and AN D = -AN A (initially, AN D = 0) [8,16], As a difference with photovoltaic contribution to the current density, drift and diffusion terms involve the full crystal charge distribution.This implies that the calculation of n for each cell requires solving a system of linear equations of size M 2 , where M= m x x m y is the number of cells in the surface.On the other hand, the time evolution of bounded charge (ANn) may be calculated locally for each cell of the crystal, and is given by Eq. ( 2 The integration procedure is made by alternatively solving n and AND distributions in time steps dt. Once the electric field of each cell is obtained, in order to calculate the evanescent field outside the material, it is necessary to use an electric field distribution along the Z-axis inside the crystal.We will assume that the electric field obtained for each cell is homogeneously distributed along the Z-axis through the whole thickness of the cell (see Fig. 1(b)).

External electric field
To calculate the external electric field in a point P, in the near vicinity of the crystal E p , each cell is now converted to filament.So a given filament with the previously obtained bounded charge ATVp is located in position (i,j) on the crystal surface, with a depth l z and a surface size {AxAy) (see Fig. 1(b)).
The contribution of this filament to E p \s given by: where n points to the charge element in the filament from the crystal surface, r 2 points to P from the surface of the filament, z 2 and pi are the projections of, normal and parallel to the surface, respectively and 9 is the angle between the direction of pi and the C-axis.The electric field in P(x,yz) is the summation of the contributions of all the Nix A^/filaments: Although the depth of the sample has no incidence in the calculation of the charge distribution, it has influence in the definition of the electric field near the surface of the crystal.Thinner samples would have sharper E-profiles near the surface than thicker ones.
Once the value of i^is obtained Eq. ( 6) directly provide the dielectrophoretic potential at P.

Results of simulations
This algorithm has been used to calculate the electric field as well as the DEP potential at any point of the surface of the crystal and near it for different representative light patterns.In the following subsections the effect of different light patterns are shown.

Sinusoidal light patterns
First, we will consider a sinusoidal light pattern parallel to the C-axis, where the light intensity varies according to the expression being k the average value of the intensity, /nthe modulation and Kthe modulus of the grating vector.This light pattern is very easy to achieve with the interference of two laser beams and the modulation can be controlled with the relative intensities of the beams.The cases of low (/n=0.1)and high (/n=0.95)modulation have been studied for different fringe periods (A).A sinusoidal grating with low modulation is a well known illumination pattern for which, analytical solutions were previously obtained (see for instance [8,16]).Figure 2 shows the internal charge distribution, ANo, and the internal electric field, E x , for a given low modulation intensity pattern along the C-axis (X direction).First, the charge distribution AAfc> and the internal electric field E, generated by an intensity grating with m = 0.1 along the C-direction are calculated.As expected for this low modulation, ANo shows an almost sinusoidal distribution with a phase shift of n/2 with the intensity pattern in accordance with previous theory.
Nevertheless, despite the clearly oscillating modulation of the electric field shown in Fig. 2, it is known [14] that, for low modulation, the electric field pattern inside the crystal should only produce an homogeneous particle trapping.Moreover recent experimental results corroborate this prediction [19].This fact can easily be explained by examining the electric field distribution and the DEP potential outside the crystal.From Fig. 3 one can notice a non-homogeneous distribution of the electric field, but a carefully study of this field shows that, although the direction of the field changes along the X-axis its module remains constant and changes only along the Z direction.This leads to a constant DEP potential along the X-axis that increases in absolute value as one moves toward the crystal surface (see Fig. 3).In fact the DEP potential qAxAyANfi 1  4ne 0 e r Fig. 2. Light intensity /with m = 0.1, internal charge tS.No an d internal electric field EQis not completely constant, it has a small curl that increases as one moves closer to the surface.This effect can be explained due to continuity considerations since the electric field must be continuous at crystal surface (inside and outside).It is well known from the standard theory of the photorefractive effect that, at high modulation, charge transport is a nonlinear process and the space charge field is not sinusoidal but only periodic with an important contribution of higher harmonics of the grating vector K, predicted by theory [16] and recently observed in Fe:LiNb03 [19].As a consequence, the space charge field localizes and becomes asymmetric [20].These predictions are in good agreement with the results obtained by our algorithm at m = 0.95.They are represented in Fig. 4, where the space charge field and the corresponding DEP potential are plotted.As expected, the charge distribution is more localized and asymmetrical and the DEP potential has well defined minima that become less pronounced as the vertical distance to the substrate increases.Its periodicity along the X-axis coincides with that of the light.About this point there has been some controversy because some authors expected a grating vector 2k.This result confirms previous analytical predictions by Villarroel et al. [14] and is in accordance with reported experiments in [15,19], Notice that the minima of DEP potential in Fig. 4(b) (i.e. the positions of trapped particles) coincide with the light intensity maxima.
As the trapping of dielectric particles is not affected by the internal electric field, in the rest of this work, only the external electric field and DEP potential will be shown.

Gaussian light intensity distribution
Avery interesting and simple case is the illumination by a single gaussian beam (of width a).
Figure 5 shows the illumination pattern and the charge distribution generated by it.It can be seen that photovoltaic current generates a migration of charge along the C-axis.The electric field and the DEP potential above the crystal surface are shown in Figs. 6 and 7, respectively.Trapping of particles will occur at the minima of the potential, which have a crescent shape, and are located symmetrically on both sides of the beam along the C-axis, slightly outer than the maxima of | AV\fc>| .As expected, the closer to the surface the better the potential shape is defined.
C axis (a) C axis (a)

Interference of two Gaussian beams
The last case studied corresponds to the interference of two gaussian beams.Actually, almost any recorded grating is made by some sort of interference of gaussian beams.The resulting intensity pattern is a sinusoidal grating, convoluted by a gaussian function of the size of the interfering beams.Figure 8 shows the results for a grating period A and m = 0.7.The spatial charge field is a mixture of the two cases seen above: charge accumulated at both sides of the illuminated region and ripped with the spatial frequency of the sinusoidal grating.Surprisingly, the central part of the pattern, where light intensity is higher and grating contrast is more evident, has very low charge variation, compared with the extremes of the light pattern.According to this spatial charge distribution, the DEP potential (Fig. 8(b)) has its maxima values at the extremes of the light pattern and, despite the high contrast of the light grating, in the center of the illumination pattern there is no significant grating of potential and no trapping is expected there.

Comparison to experiment
In this section the results obtained will be compared with some representative experiments of photovoltaic micro-particle trapping.Experiments have been carried out in 1 mm thick x-cut congruent LiNbC>3 crystal doped with iron (0.1% wt) in order to have a strong photovoltaic effect.Single beam and holographic configurations have been tested.After illumination with 532 nm laser light, chalk (mainly CaCOs) micrometric particles, with a size range of 1 -lOjiim, were either sprayed from beneath the bottom surface of the sample or immersed in a hexane solution where the sample was later introduced.For more details on the experimental method see [15,19].Then, the particle distribution attached to the surface was visualized either by an optical photograph or a micro-photograph.Figurem9 shows the images of particle patterns obtained with three different light distributions, corresponding to the three cases analyzed in section 4. Figure 9(a) has been obtained by particle deposition after illumination with a single incoherent beam having an approximate gaussian light profile.Figures 9(b) and 9(c) correspond to light patterns generated by the interference of two approximate plane waves and of two gaussian beams respectively.For both interferograms, the two light beams had the same intensity so that the modulation factor, m, is close to 1.The fringe period is 30 fim in pattern of Fig. 9(b) and lOOjiim in pattern of Fig. 9(c).It is clearly appreciated that the general morphology of the experimental patterns is in very good accordance with our predictions for the DEP potential of the model.In Fig. 9(a) particles are trapped at the boundary of the light pattern at both sides of the C-axis (horizontal axis in the Fig. 9) following the same distribution found in Fig. 7.For a sinusoidal light pattern, we obtain a periodic trapping with the same period than lighting, as predicted by our model (see Fig. 4(b)).Finally, in Fig. 9(c) the particle pattern is also periodic with the same periodicity as lighting, but the central region is free of particles in accordance with a DEP potential nearly free of local minima (see Fig. 8(b)).Conversely, the particles are trapped at both sides of the region with pronounced minima in the DEP potential.
Also, as pointed out in section 4.1.1,for a low modulation sinusoidal light pattern, soft homogeneous trapping is observed.
As a final step, we have compared our numerical results with a more general two dimensional experiment taken from [15] as reported in [19].In that work, the sample was illuminated with a light pattern of circular fringes obtained from the interference between plane and spherical waves (Fresnel-type pattern).The pattern diameter was 9 mm-long and the ring separation between 250 and 500 /im (Fig. 10(a)).After the pattern was created, chalk particles where sprayed from beneath the exposed surface creating the chalk pattern on the surface of the crystal of Fig. 10(c).Figure 10(b) shows the DEP potential obtained at different distances from the crystal surface.Again, very good agreement between theory and experiment is obtained.The larger trapping on the outer part of the fringes, probably comes from the attraction of the particles at larger distances from the sample.There, only two crescent shape minima appear in the DEP potential, at both sides of the illuminated zone.It can also be observed how, despite the illumination of the rings being homogeneous, both, the DEP potential and the chalk pattern show no trapping in the central line perpendicular to the C-axis.

Conclusions
An algorithm to calculate the space charge field and the electric field generated by an arbitrary illumination pattern on the surface of a X-cut photorefractive crystal has been developed.This electric field has been used to calculate the DEP potential generated by illumination on the surface of the crystal and nearby.This algorithm has been tested using the physical parameters given for iron doped lithium niobate for several, common and exotic, illumination patterns.The obtained results have been compared with existing experiments with a very good agreement.So, this method offers a way to predict the trapping capability of an illumination pattern and can be used to determine the best one to get a desired particle trapping design.

Fig. 1 .
Fig. 1.(a) Two-dimensional representation of the crystal.Y-direction along C-axis.(b) Reference system and charge filaments for external field Ef, calculation.

Fig. 8 .
Fig. 8. Gaussian beam interference with grating period A and m=0.7:(a) intensity (I) and ANJJ; (b) dielectrophoretic potential in air at several distances from the crystal surface.

Fig. 10 .
Fig. 10.Trapping of Fresnel-type pattern: (a) light intensity, (b) DEP potential over the surface of the crystal and (c) experimental trapping of CaCC>3 particles. .

Table 1 .
Photovoltaic parameters used in calculations