Computer generation of optimal holograms for optical trap arrays

: We propose a new iterative algorithm for obtaining optimal holograms targeted to the generation of arbitrary three dimensional structures of optical traps. The algorithm basic idea and performance are discussed in conjunction to other available algorithms. We show that all algorithms lead to a phase distribution maximizing a speciﬁc performance quantiﬁer, expressed as a function of the trap intensities. In this scheme we go a step further by introducing a new quantiﬁer and the associated algorithm leading to unprecedented efﬁciency and uniformity in trap light distributions. The algorithms performances are investigated both numerically and experimentally.


Introduction
A mesoscopic object can be stably trapped in three dimensions by a tightly focusing single laser beam [1].Computer-generated holograms (CGH) displayed on liquid crystal spatial light modulators (SLM) offer a convenient way of producing large three dimensional arrays of optical traps [2,3,4,5,6].The ability to dynamically manipulate matter at the meso-scale opened the way to a wide range of applications in the physical and biological sciences [7].
In the usual implementation of holographic optical trapping an SLM is used to modulate the complex amplitude of a laser beam in the back focal plane of an optical system terminating with a high magnification microscope objective.On each pixel of the SLM, the incident beam undergoes both a phase and polarization modulation whose entity can be usually controlled with an eight bit depth.SLMs, being usually based on liquid crystals displays, do not allow an independent modulation of both phase and polarization.On the contrary, for each pair of incident and outgoing polarization directions, only a one dimensional path on the complex modulation plane can be accessed by varying the value of the displayed gray level.In optical tweezers applications it is usually more convenient to use all the available power and impose a phase only modulation.Unfortunately the task of searching for a phase only hologram producing the desired intensity distribution is not straightforward.To this aim many algorithms have been proposed in literature, each one having advantages and drawbacks.In particular, as discussed in the following section, when dealing with highly symmetric trap geometries one has to abandon at least one of the following requirements: efficiency, uniformity, symmetry.In this paper we introduce an iterative procedure which allows to obtain arbitrary three dimensional traps arrangements with a theoretical efficiency typically greater than 90% and arbitrarily high uniformity.We derive our algorithm as a further step in the sequence of previously available algorithm, each one interpreted as a maximization task of some performance quantifier written in terms of trap intensities.The superior performance of the algorithm is demonstrated both numerically and experimentally for the case of a 10 × 10 square grid of point traps.

Algorithm description and performance
We will now briefly review the most commonly used algorithms for CGH, limiting ourselves to optical tweezers applications where target intensities are arbitrary three dimensional collections of high light intensity points (point traps).We will assume a uniform illumination on the SLM and call u j = |u| exp(iφ j ) the complex amplitude of electric field reflecting off the jth pixel, where φ j is the corresponding phase shift.The total energy flux through the SLM is given by W 0 = cε 0 N|u| 2 d 2 /2 where N is the total number of pixels and d 2 is the pixel's surface area.We can use scalar diffraction theory to propagate the electric field complex amplitude from the jth pixel surface to the location of the mth trap in image space [9].Summing up the contributions from all the N pixels we obtain the complex amplitude v m of electric field on trap m: where x j ,y j are the jth pixel's coordinates on the back focal plane (SLM) and x m , y m , z m are the mth trap coordinates referred to the front focal plane (Fig. 1).We can easily generalize the ∆ m j to add orbital angular momentum to trapping beams [5].To make the notation more compact we introduce the adimensional quantity V m : whose physical meaning can be understood noting that I m = |V m | 2 measures the energy flux in units of W 0 flowing through an area f 2 λ 2 /(Nd) (the area of a diffraction peak) centered at the mth trap site.For z m = 0 V m , corresponds to the discrete Fourier transform of e iφ j evaluated at the spatial frequencies (x m /λ f , y m /λ f ).
Our task here is to search, for a given set of ∆ m j , for the best choice of φ j s to maximize the modulus of V m on all traps.We will use as a benchmark to compare different strategies, the task of computing an N = 768 × 768, 8 bit hologram aiming at a target intensity of M = 100 traps arranged on a 10 × 10 square lattice located in the Fourier (z m = 0) plane.The performance of different strategies is quantified by three parameters: efficiency (e), uniformity (u) and percent standard deviation (σ ) where .. denotes the average over trap index m.Let's start from the trivial case of one single trap, M = 1.The choice here is easily found by setting φ j = ∆ 1 j which makes all terms in the sum (3) real and equal to 1/N, thus giving |V 1 | 2 = 1.When M > 1 we have to seek for a compromise between the M different choices φ j = ∆ m j (one choice for each m value) that would divert all energy on trap m.One of the fastest routes is the random mask encoding technique (RM) [17].The compromise is obtained by setting: where m j is a number between 1 and M randomly chosen for each j.The technique is very fast, and performs remarkably good as far as uniformity is concerned.However the overall efficiency can be very low when M is large.In fact, on average, for each m only N/M pixels will interfere constructively, all the others giving a vanishing contribution.Therefore |V m | 2 ≃ 1/M 2 and e ≃ 1/M which can be significantly smaller than one when M is large.In the present case, where M = 100, we numerically obtained u = 0.58 but e = 0.01 = 1/M.
Another commonly used compromise is obtained as follows.We cannot make simultaneously all V m real and equal to 1/ √ M as in the M = 1 case but we can try to maximize the real part of ∑ m V m with respect to φ j .The stationary points are easily obtained imposing the condition of a vanishing gradient: whose solutions are given by: For the stationary point to be a local maximum the corresponding Hessian matrix has to be negative defined.In the present case the Hessian is purely diagonal, and when evaluated on the stationary points reads: The maximum condition is obtained when all the n j are set to 0 and therefore: which can be read as the phase of the sum of single trap holograms and it's usually called "superposition of prisms and lenses" (S) [2,3].The S algorithm, though slower than RM (due to the extra N arg function evaluations), gives order one efficiencies but very poor uniformities.In this case, infact, though efficiency raises to e = 0.29 uniformity is only u = 0.01.Moreover, when, as for the square lattice, highly symmetrical trap geometries are sought, a consistent part of energy is diverted to unwanted ghost traps.A better compromise is obtained if we only try to maximize the sum of the amplitudes of V m projected on randomly chosen directions in complex plane.In other words we seek a maximum of ∑ m Re {V m exp(−iθ m )} where θ m are random numbers uniformly distributed in [0, 2π].In this case we obtain: (10) which is the phase of the linear superposition of single trap holograms with coefficients of unit modulus and random phase.This last choice, usually called Random Superposition (SR) [10], is of the same computational cost as S, but results in a much better performance than S, though trap intensities still vary a lot (e = 0.69, u = 0.01).
We want to stress here that, when dealing with low symmetry geometries, SR holograms can also produce good uniformity levels and no further refinement is needed.If precise trap positioning is not an issue one can deliberately reduce the pattern symmetry by adding a small amount of random displacement to trap locations as demonstrated in [11].
To go a step further from SR we have to release the constraint given by the randomly chosen phases θ m and try to maximize ∑ m |V m | allowing for any possible value for θ m .Again differen- tiating with respect to φ j we obtain the stationary points: This time the Hessian computed on the stationary point is not purely diagonal: However the non diagonal terms are 1/N smaller than the diagonal ones.It can be shown that such a perturbation will only affect the sign of one eigenvalue at most [12].When N is very large we can neglect this eventuality and call the stationary point a maximum.In this case φ j are obtained as the phase of the linear superposition of single trap holograms with coefficients of unit modulus and a phase given by the phase of V m , that is the field produced by the φ j themselves on trap site m.It's now impossible to have the φ j in an explicit form given the implicit dependence of V m on φ j .A possible route is that of starting with a guess for φ j , i.e. the one obtained from SR, and use (14) in an iterative procedure.This algorithm is called Gerchberg-Saxton (GS) [13,14] and it converges after a few tens of iterations.In particular after 30 iterations we obtained e = 0.94, and u = 0.60.The reason why the algorithms so far discussed usually result in poor uniformities can be understood if one observes that we are always concerned with the maximization of the sum of the amplitudes of V m , having no bias towards uniformity.Such a bias is present when we seek for a maximum in a quantity like ∏ m |V m | or equivalently ∑ m log |V m |.By differentiating the biased function with respect to φ j we obtain: It is easy to show that also in this case the Hessian matrix is diagonal (in the limit of large N) and negative defined at the stationary point: If we seek the solution of ( 16) by an iterative procedure we obtain the Generalized Adaptive Additive algorithm (GAA) [5].With the choice ξ = 0.5 GAA produces a uniformity improvement u = 0.79 with the same efficiency e = 0.93 as GS.
The last algorithm we review achieves the maximization of a gain function by a direct search.Starting from a good guess such as one obtained from SR we pick up one pixel at random and cycle through all the P = 256 gray levels looking for an improvement (increase) in the gain function: This is the Direct Search algorithm (DS) [15,16].As suggested in [16] when starting from a SR hologram and setting f = 0.5, the algorithm achieves a perfect uniformity (u = 1.00) after 1.3 × N steps of computational cost scaling as M × P, though the overall efficiency is diminished to e = 0.68.Better holograms can be obtained by giving more bias to efficiency ( f =0.25) and waiting for a substantially longer time (∼ 10 N steps that is about a hundred times longer than GS).However, we observed that reducing the number of gray-levels P to just 8 can significatively reduce (by a factor 32) the computational cost without affecting performance too much (see [15] for a systematic exploration of parameter space).With the previous parameter choice and with 8 grey levels we obtained e = 0.84 and u = 1.00 after 7 N steps that is still about 3 times longer than GS.A this point the whole hologram has been reduced to 3 bit and a comparison with other algorithms working with full 8 bit depth is out of purpose.So far, we reviewed the most common strategies for generating phase only holograms for optical tweezers applications.Each one of the presented methods is not fully satisfying regarding at least one parameter between efficiency and uniformity.
We will now go a step further introducing a new iterative GS like algorithm having optimal performances with respect to both efficiency and uniformity.We introduce the M extra degrees of freedom w m , that maximize the weighted sum ∑ m w m |V m | with the constrain that |V m | are all equal.By differentiating with respect to φ j , we obtain the maximum condition: Again, the above formula expresses φ j in an implicit form, this time containing also the un- known weights w m .Starting from a SR guess for φ j and setting w m = 1, the iteration proceeds as follows: In other words, at each step we adjust the weight w m in such a way to reduce |V m | deviations from the average |V | .The above procedure converges, with a speed typical of GS and GAA, to a hologram having the almost optimal performance e = 0.93, u = 0.99.We will refer to this new algorithm as weighted Gerchberg-Saxton or GSW.The optimization progress for GSW is reported in detail in Fig. 2 and shows how we can efficiently use it to obtain above 90% efficiency and uniformity in only 10 steps.The performance of our algorithm remains the highest when we turn to three dimensional lattices.For example a 3 × 3 × 3 simple cubic lattice can be generated with an efficiency of e = 0.91 and a u = 0.99 uniformity.The idea of aiming at a slightly modified target intensity distribution in GS optimizations has been also proposed in a different algorithm targeted to two dimensional beam shaping tasks [18].
Table 1 summarizes the results of the benchmark test described above, clearly showing the superior performance of GSW with respect to the currently available algorithms.

Experimental results
The experimental setup is schematically represented in Fig. 3.A TEM00 mode beam from a diode pumped, λ =532 nm, 2 W laser (LaserQuantum Opus) is expanded and reflected off a liquid crystal (45 • twisted nematic) Spatial Light Modulator (Holoeye LCR-2500).The SLM is placed between two linear polarizers adjusted to obtain a phase mostly modulation (P = −4 • , A = 20 • , where angles are measured from the vertical direction and rotating clockwise when viewing from the beam direction).The phase modulated wavefront is imaged onto the exit pupil of a 100x NA 1.4 objective lens mounted in a Nikon TE2000-U inverted optical microscope.The light distribution in the objective focal plane is deduced from the brightness of the laser light diffusing water-coverglass interface imaged on a CCD (Prosilica).The spacing between traps in the target array has been set to 3 µm.We obtain experimental values for I m by summing the values of pixels inside a 1µm × 1µm square area centered on the peak center of mass.Such a determination of I m is affected by many sources of uncertainty: i) the coverglass plane might not be exactly parallel to the traps plane, ii) the diffusing power of the interface might be non uniform, ii) the imaging with coherent and polarized light might be non uniform.To correct for these uncertainties, we measured an ensemble of I m obtained from high uniformity algorithms (GSW, DS) and calculated starting from independent SR holograms.For each m, we thus obtained a normalization factor as the average of the measured I m s.The corrected I m for the algorithms described in section 2 are reported, together with the raw image data, in Fig. 4.
We excluded RM algorithm in this comparison due to its bad performance in this case.All algorithms perform as calculated theoretically, in particular GSW shows a greater efficiency and uniformity.It is clear from Fig. 4 that the low efficiency of DS is due to light going into ghost traps located in the nearby sites of the square lattice.To better compare with theory we report the parameters e,u,σ in Table 2. Efficiencies have been multiplied by a constant factor chosen to give the theoretically expected value for GSW.The agreement with the expected performance is very good, but we never observe a perfect uniformity (u ≃ 1), probably for some residual error in the determination of I m .

Conclusion
In conclusion, we reviewed most of the available algorithms for the generation of phase only holograms targeted to optical tweezers applications.With the exception of RM, we showed that all algorithms lead to a phase distribution representing the maximum point of some performance quantifier expressed as a function of the trap intensities I m .In this scheme, we went a step further by introducing a new algorithm leading to unprecedented efficiency and uniformity in trap light distributions.The algorithm converges in a few tens of iterations of N × M computational cost.We compared, both numerically and experimentally, the performances of investigated algorithms in producing a 10×10 square grid target.The obtained results demonstrate that the proposed algorithm allows to achieve almost perfect efficiencies and uniformities using phase only holograms and a modest computational time.

Fig. 1 .
Fig. 1.Schematic representation of Fourier optics propagation from SLM plane (back focal plane) to the front focal plane of the optical system.

Fig. 4 .
Fig. 4. Experimental determination of light distribution on trapping plane.For each algorithm we report the raw image of laser spots (left) together with the corrected (see text) values for trap intensities (I m ) represented in the gray level of square tiles centered at the corresponding lattice site.

Table 1 .
Summary of theoretical performances of investigated algorithms.The target trap structure is a 10 × 10 square grid.Column 2 contains a 100 × 100 detail of the total 768 × 768 hologram.Performance parameters after K (column 6) iterations are reported in columns 3,4,5.Computational cost scaling is reported in column 6 where: M=number of traps, N=number of pixels in hologram, K=number of iterations, P=number of gray levels (256 here).

Table 2 .
Summary of experimental performances of investigated algorithms.