Minimizing intensity fluctuations in dynamic holographic optical tweezers by restricted phase change

We present a method for reducing intensity fluctuations that typically occur when a spatial light modulator is updated between consecutive computer generated holograms. The method is applicable to most iterative hologram generating algorithms and minimizes the average phase difference between consecutive holograms. Applications with high stability requirements, such as optical force measurement with holographic optical tweezers, should benefit from this improvement. ©2010 Optical Society of America OCIS codes: (090.1760) Computer holography; (120.5060) Phase modulation; (140.7010) Laser trapping; (230.6120) Spatial light modulators; (350.4855) Optical tweezers or optical manipulation; (090.2890) Holographic optical element; (170.4520) Optical confinement and manipulation References and links 1. A. van der Horst, and N. R. Forde, “Calibration of dynamic holographic optical tweezers for force measurements on biomaterials,” Opt. Express 16(25), 20987–21003 (2008). 2. C. O. Mejean, A. W. Schaefer, E. A. Millman, P. Forscher, and E. R. Dufresne, “Multiplexed force measurements on live cells with holographic optical tweezers,” Opt. Express 17(8), 6209–6217 (2009). 3. A. Farré, A. van der Horst, G. A. Blab, B. P. B. Downing, and N. R. Forde, “Stretching single DNA molecules to demonstrate high-force capabilities of holographic optical tweezers,” J. Biophoton. 3(4), 224–233 (2010), doi:10.1002/jbio.200900107. 4. D. Preece, R. Bowman, A. Linnenberger, G. Gibson, S. Serati, and M. Padgett, “Increasing trap stiffness with position clamping in holographic optical tweezers,” Opt. Express 17(25), 22718–22725 (2009). 5. M. Reicherter, S. Zwick, T. Haist, C. Kohler, H. Tiziani, and W. Osten, “Fast digital hologram generation and adaptive force measurement in liquid-crystal-display-based holographic tweezers,” Appl. Opt. 45(5), 888–896 (2006). 6. 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(9), 608–610 (1999). 7. E. R. Dufresne, G. C. Spalding, M. T. Dearing, S. A. Sheets, and D. G. Grier, “Computer-generated holographic optical tweezer arrays,” Rev. Sci. Instrum. 72(3), 1810–1816 (2001). 8. D. G. Grier, “A revolution in optical manipulation,” Nature 424(6950), 810–816 (2003). 9. L. P. Ghislain, N. A. Switz, and W. W. Webb, “Measurement of piconewton forces using a simple optical force microscope,” Biophys. J. 66, A278 (1994). 10. J. Liesener, M. Reicherter, T. Haist, and H. J. Tiziani, “Multi-functional optical tweezers using computergenerated holograms,” Opt. Commun. 185(1-3), 77–82 (2000). 11. T. Haist, M. Schonleber, and H. J. Tiziani, “Computer-generated holograms from 3D-objects written on twistednematic liquid crystal displays,” Opt. Commun. 140(4-6), 299–308 (1997). 12. D. Engström, A. Frank, J. Backsten, M. Goksör, and J. Bengtsson, “Grid-free 3D multiple spot generation with an efficient single-plane FFT-based algorithm,” Opt. Express 17(12), 9989–10000 (2009). 13. J. E. Curtis, B. A. Koss, and D. G. Grier, “Dynamic holographic optical tweezers,” Opt. Commun. 207(1-6), 169–175 (2002). 14. M. W. Farn, “New iterative algorithm for the design of phase-only gratings,” in Proc. SPIE(1991), pp. 34–42. 15. R. Di Leonardo, F. Ianni, and G. Ruocco, “Computer generation of optimal holograms for optical trap arrays,” Opt. Express 15(4), 1913–1922 (2007). #126429 $15.00 USD Received 5 Apr 2010; revised 7 May 2010; accepted 7 May 2010; published 12 May 2010 (C) 2010 OSA 24 May 2010 / Vol. 18, No. 11 / OPTICS EXPRESS 11250 16. E. Eriksson, S. Keen, J. Leach, M. Goksör, and M. J. Padgett, “The effect of external forces on discrete motion within holographic optical tweezers,” Opt. Express 15(26), 18268–18274 (2007). 17. E. Marom, and N. Konforti, “Dynamic optical interconnections,” Opt. Lett. 12(7), 539–541 (1987). 18. M. Johansson, S. Hård, B. Robertson, I. Manolis, T. Wilkinson, and W. Crossland, “Adaptive beam steering implemented in a ferroelectric liquid-crystal spatial-light-modulator free-space, fiber-optic switch,” Appl. Opt. 41(23), 4904–4911 (2002). 19. E. Hällstig, J. Öhgren, L. Allard, L. Sjöqvist, D. Engström, S. Hård, D. Ågren, S. Junique, Q. Wang, and B. Noharet, “Retrocommunication utilizing electroabsorption modulators and nonmechanical beam steering,” Opt. Eng. 44(4), 045001 (2005). 20. R. W. Gerchber, and W. O. Saxton, “Practical algorithm for determination of phase from image and diffraction plane pictures,” Optik (Stuttg.) 35, 237 (1972). 21. F. Wyrowski, and O. Bryngdahl, “Iterative fourier-transform algorithm applied to computer holography,” J. Opt. Soc. Am. A 5(7), 1058–1065 (1988). 22. D. Engström, G. Milewski, J. Bengtsson, and S. Galt, “Diffraction-based determination of the phase modulation for general spatial light modulators,” Appl. Opt. 45(28), 7195–7204 (2006). 23. E. Hällstig, J. Stigwall, T. Martin, L. Sjöqvist, and M. Lindgren, “Fringing fields in a liquid crystal spatial light modulator for beam steering,” J. Mod. Opt. 51(8), 1233–1247 (2004). 24. J. L. Harriman, A. Linnenberger, and S. A. Serati, “Improving spatial light modulator performance through phase compensation,” Proceedings of the SPIE The International Society for Optical Engineering 5553, 58–67 (2004). 25. X. Xun, and R. W. Cohn, “Phase calibration of spatially nonuniform spatial light modulators,” Appl. Opt. 43(35), 6400–6406 (2004). 26. J. Otón, P. Ambs, M. S. Millán, and E. Pérez-Cabré, “Multipoint phase calibration for improved compensation of inherent wavefront distortion in parallel aligned liquid crystal on silicon displays,” Appl. Opt. 46(23), 5667–5679 (2007). 27. A. Jesacher, A. Schwaighofer, S. Fürhapter, C. Maurer, S. Bernet, and M. Ritsch-Marte, “Wavefront correction of spatial light modulators using an optical vortex image,” Opt. Express 15(9), 5801–5808 (2007). 28. J. E. Curtis, C. H. J. Schmitz, and J. P. Spatz, “Symmetry dependence of holograms for optical trapping,” Opt. Lett. 30(16), 2086–2088 (2005).


Introduction
Optical force measurement (OFM) using holographic optical tweezers (HOT) is a promising and rapidly emerging technique providing a powerful and versatile platform for novel research in the life sciences [1][2][3][4][5].It combines the ability of simultaneous independent manipulation of a large number of objects [6][7][8] with the possibility of precise measurements of forces on a biologically relevant scale [9].In a HOT setup, spot positioning is accomplished by diffraction of the trapping beam by a spatial light modulator (SLM).Dynamic beam steering is realized by displaying a sequence of phase patterns on the SLM, each yielding a desired spot configuration.Thus, each phase pattern, here referred to as a computer generated hologram (CGH), in the sequence needs to be calculated in such way that it produces the desired spot pattern.From the task the HOT is to perform, a set of requirements such as diffraction efficiency, uniform intensity among the spots, and allowed calculation time are given.Depending on these requirements the CGHs can be calculated with various methods.The fastest method is to use the resulting phase after complex addition of a set of simple CGHs, each corresponding to one of the desired spots [10].Such CGHs generally suffers from non-uniform intensity distribution among the spots and additional ghost spots at undesired positions.CGHs of higher quality are obtained with iterative algorithms based on Fraunhofer propagation [11,12] or on Fresnel propagation [13][14][15].
Examples of basic tasks for a HOT system using multiple beams are: to turn off a single spot in order to release a trapped object; add a new spot in order to trap an additional object; and to move one or more spots in small steps.Despite that the changes in the spot patterns may be relatively small in the three tasks, strong intensity fluctuations generally occur during the update of the SLM.An example of this is shown in Fig. 1, where one of 25 spots has been removed in the second CGH.The cause of the fluctuations is that during the update from one CGH to another, the SLM will undergo a sequence of more or less uncontrolled states until all pixels have reached their new settings.This typically causes light to be temporarily directed away from the desired spots.Such intensity fluctuations result in corresponding fluctuations in the trap stiffness.As illustrated in Fig. 1(d), the intensity may even vanish in some of the spots during the SLM update, with the risk of losing one or more of the trapped objects as a consequence.If the SLM is continuously updated to a new CGH at a time shorter than its response time, the SLM pixels may even remain in non-favorable states (i.e., something between the consecutive CGHs) so that the spot intensities will never reach their intended values.When combined with OFM, the fluctuations are clearly a problem as very high stability and predictability of trap stiffness is required.Said intensity fluctuations have previously been studied for the case where a single spot was moved using CGHs consisting of blazed gratings [16].For this case it was shown that the magnitude of the intensity loss is related to the average difference in phase setting of the SLM pixels between consecutive CGHs.The fluctuations could therefore only be minimized by using a very small step size for the moving spot which hence increased the similarity of the consecutive CGHs.
In this work we have studied this effect for more complex spot patterns consisting of three or more spots.We present a method in which the maximum allowed phase change between two consecutive CGHs is constrained for each pixel.The method, which we call "Restricted Phase Change" (RPC), forces the algorithm to converge to a CGH with much stronger correlation to the previous one than with conventional methods.With this method, the amplitude of the intensity fluctuations can be made negligibly small for all three basic HOT tasks mentioned above.When moving spots, there is no need to reduce the step size of the moving spots.Thus, faster movements can be allowed.
Obvious advantages of using SLMs for beam steering are that individual 3D positioning of the spots, abrupt changes of the spot pattern, individual spot shaping, and combinations thereof are readily accomplished.Due to these benefits, SLMs are used for non-mechanical beam steering not only in HOT systems but also in applications such as optical interconnects [17,18] and free-space communication systems [19].Although this work focuses on the use of SLMs for HOT applications, the issues tackled are therefore of importance in any SLM-based beam steering system using multiple beams.

Theory
Before describing our method in Section 2.2, a brief description of a general GS algorithm will be given.

Hologram generating algorithms
Some of the most popular algorithms for generating optimized phase only holograms are variations of the Gerchberg-Saxton algorithm (GS) that was first introduced in 1972.It was originally used to retrieve the phase distribution of some incident field by observing its intensity distribution and the intensity distribution of its corresponding diffraction pattern [20].For hologram generation, it is instead used to determine the incident phase distribution necessary to obtain some desired intensity distribution in the output field.The incident field usually has a fixed amplitude distribution determined by the beam profile of the laser while the phase may be set arbitrarily, only restricted by the pixelation and bit depth of the SLM.
A large number of algorithms that produce spot generating phase holograms can be classified as modified GS algorithms [7,13,14,21].Common principles for such algorithms are presented in Fig. 2: A complex incident field is created by combining the laser amplitude profile, ( , )  where max I and min I refers to the intensity of the brightest and the weakest spot, respectively.GS algorithms are well suited for execution on graphics processing units (GPU) as numerical propagation methods generally are highly parallelizable.In our implementation, highly optimized 512 × 512 CGHs for 3D positioning of 20 spots are generated in 10 ms using a consumer grade GPU (NVIDIA GeForce GTX 285).

Reducing intensity fluctuations during SLM update
GS algorithms usually result in CGHs with high diffraction efficiency and virtually perfect uniformity.However, as previously mentioned, unwanted fluctuations of the spot intensities may occur during the update of the SLM.In [16], Eriksson et al. showed that the magnitude of the fluctuations are closely related to the average pixel phase difference between consecutive CGHs.As part of this work we have investigated how this entity is influenced by the choice of 0 ( , ) x y φ .We have compared using a random initial guess, here referred to as ( , ) Rand x y φ , for the GS algorithm and using the phase distribution of the previous CGH, here referred to as 1 ( , ) CGH x y φ .For convenience, we will refer to the two methods as GS-Rand and GS-CGH1.To force the algorithm to converge to a CHG with lower average pixel phase difference, we also introduce a novel method that limits the maximum phase difference allowed for each pixel between consecutive CGHs.This is done by introducing an additional step in the GS iteration cycle as shown in Fig. 3: After the propagation back to the plane of the SLM, the where α is a constant set by the user which determines the maximum allowed phase change for each pixel.Thus, we call our method RPC (restricted phase change) and the implementation used here is referred to as GS-RPC.The choice of α indirectly determines the number of pixels that remains unaltered in the new CGH.If α is set to 1, the method is equivalent to GS-CGH1, while setting α to 0 will result in a CGH identical to the previous one.The optimal choice of α may vary depending on the application.Guidelines for selecting an appropriate α will be given in section 4.2.
The technique was implemented in two versions of the GS algorithm, one that used fast Fourier transforms (FFTs) for propagation between the plane of the SLM and the spot positions and one that used Fresnel propagation.In HOT applications, trap positions are constrained to an equidistant grid in a single plane using the former version while the latter allows for arbitrary 3D positioning of traps.

Experimental methods
The intensity fluctuations were measured experimentally by imaging the far field of a spatial light modulator directly onto the sensor chip of a CMOS camera in a free space setup and by observing the position of trapped particles in a constant fluid flow in an optical tweezers setup.The CGHs used in all experiments were generated with a weighted FFT-based version of the Gerchberg-Saxton algorithm.
In the following sections, distances in the trapping plane will be given normalized to the size of a diffraction limited spot in this plane.The spot size is here defined as fλ/D, where f is the focal length of a Fourier transforming lens, λ is the wavelength and D is the beam width in the back focal plane of the Fourier transforming lens.The unit "spot size" is thereby independent of the optical system used and corresponds to the lattice spacing given by the FFT.

Free space experiments (BNS SLM)
As shown in Fig. 4, a 125 mm Fourier transforming lens was used to image the far field of the SLM (BNS 512 × 512 SLM System, electrically addressed) directly onto the chip of a CMOS camera (Photonfocus MV-D1024E-160-CL).The SLM was illuminated by a solid state laser (Laser Quantum Ventus VIS 532nm <500mW) that was expanded to match the size of the phase pattern on the SLM.To make sure the SLM frames were encoded correctly, a diffraction based characterization method [22] was used to obtain the relation between the externally controlled setting of a pixel and the resulting (stationary) value for the phase modulation.To reduce undesired effects caused by the peripheral regions of the slightly curved SLM backplane, only the central 256 × 256 pixels were used in the experiments.The remaining pixels were set to a binary grating deflecting the light to positions outside the measurement area.Images of the diffraction spots were transferred to a computer at a rate of 200 frames per second.The power contained in an area of 15 × 15 camera pixels centered at each spot was calculated and used as the value for the power in the spot in the analysis.The used CGHs formed 9 spots located among 10 equally separated positions on a circle with a radius of 70 spot sizes centered on the optical axis.In each experiment, the SLM was updated to a CGH with the missing spot at a different position.CGHs were generated using each of the methods, GS-Rand, GS-CGH1, and GS-RPC.Results are presented in Fig. 6.

Optical tweezers experiments (Hamamatsu SLM)
The optical tweezers setup, illustrated if Fig. 5, was built around a motorized inverted epifluorescence microscope (DMI6000B, Leica Microsystems).The laser beam (1070 nm, IPG Photonics) was first magnified by an afocal telescope (f1 = 50 mm, f2 = 200 mm) to match the beam diameter to the full width of the SLM (X8267-15, Hamamatsu Photonics, optically addressed).A second afocal telescope (f3 = 375 mm, f4 = 100 mm) was used to image the SLM plane onto the back focal plane of the microscope objective (100 × , NA 1.3).The magnification of the second telescope was chosen so that the output beam slightly overfilled the back aperture of the microscope objective.Bright field images were taken at a rate of 500 frames per second using a CMOS camera (Prosilica EC1280) and saved for analysis.Position data were obtained by applying a center of gravity algorithm to a small area of the images containing the trapped object after the background signal had been eliminated by thresholding.
Using each of the methods, GS-Rand, GS-CGH1 and GS-RPC, series of 10 CGHs were generated.Spots were positioned on a 5 × 5 grid, with 16 spot sizes of separation, centered on the optical axis of the system.In every second CGH, one spot, at varying positions, was removed from the grid and then replaced.For GS-RPC, α = 0.86 was used.A single 2.56µm silica bead (Bangs Laboratories, Inc.) was trapped in the lower left spot and an external force was applied by moving the microscope stage at a constant speed.The position of the trapped bead was measured while the sequences of 10 CGHs were sent to the SLM at a rate of 1 frame per second.The experiments were repeated with different speeds of the microscope stage.The trapping laser was set to 1W of output power.Results are shown in Fig. 7.

Free space experiments (Hamamatsu SLM)
The spot patterns used in the experiments described in section 3.2 were also observed with a CMOS camera (Photonfocus MV-D1024E-160-CL) in a free space configuration.A mirror was then inserted after the first lens in the second telescope, L3 Fig. 5, to project the far field of the SLM directly onto the camera chip.Images were captured at a rate of 250 frames per second while the CGH sequences were displayed on the SLM at a rate of 1 frame per second.Results are shown in Fig. 1 and Fig. 8.

Simulations
The intensity of spots during the SLM update was calculated using Fresnel propagation of holograms simulating intermediate states of the SLM.These intermediate holograms, ( , , ) x y t φ , must be generated under some assumption about the dynamics of the SLM pixels.In a first order approximation that later was evaluated, we assumed that pixel values change at a constant rate, linearly proportional to the difference between their initial and final phase values: ( )

CGH CGH CGH t x y t x y x y x y t T T
Here, T is the total duration of the update, ( , , ) In order to investigate the effect of the phase limiting parameter, α, on the intensity fluctuations during update as well as the quality of the resulting holograms, a large number of simulations were made for three different spot arrangements, each with three different step lengths for moving spots: 1. 25 spots arranged on a square grid centered on the optical axis with a separation of 10 times the spot size.The lower left spot was moved a total distance of 100 spot sizes in negative x-direction in steps of 1, 0.5 and 0.25 spot sizes.
2. 4 spots arranged and moved as above in steps of 1, 0.5 and 0.25 spot sizes.
3. 4 spots randomly distributed, one spot was moved 1000 steps in random directions in steps of 1, 0.5 and 0.25 spot sizes.
Intensities of the stationary, non-moving, spots were simulated for 20 intermediate steps for each update.The minimum and maximum spot intensity during each update were calculated as well as the diffraction efficiency of each hologram after update.The diffraction efficiency was calculated by summing the optical power in the desired spots and dividing by the total optical power in the SLM plane.The simulations were done for 51 equally distributed values of α ranging from 0.5 to 1.Each set of simulations was repeated 100 times with different initial phase distributions and, for case 3, different initial spot positions.The average and standard deviation of said variables were calculated for each value of α and for each scenario.Each sequence was preceded by the generation of a CGH with no restriction (α = 1) using the initial spot positions.The CGHs used in these simulations a size of 512 × 512 pixels and were generated with a weighted Fresnel propagation based version of Gerchberg-Saxton as described in [15], modified with RPC and using 20 iterations.Results from these simulations are presented in Fig. 9.
In order to verify the validity of our assumption regarding the dynamics of the SLM pixels during update, simulations were also made for the same set of CGHs as used in the experiments presented in Section 3.1.Results are presented in Fig. 6.  dynamics of the SLM-pixels during update appears to be slightly different from our first approximation.In order to more accurately mimic the dynamics of the SLM, various simulations were made where more complex pixel dynamics were assumed.Although the shape of some of those simulations more closely resembled experimental data, the resulting values of the parameters of interest, such as the maximum intensity fluctuation, were no different from those obtained using our original simulations which, following the principle of minimum assumptions, are presented here.

Measured and simulated temporal variation of spot power
As seen in Fig. 6, the theoretical uniformity of the spot intensities before and after the SLM update is very high while the measurements yield a less uniform distribution.This is mainly due to imperfections of the used SLM.Two common examples of imperfections that may cause a slightly distorted phase pattern, and consequentially decreased uniformity, are i) damping of high spatial frequencies in the phase pattern due to physical limitations of the liquid crystal material used to obtain the phase modulation and ii) a phase response that varies over the SLM surface.To some extent these imperfections can be compensated for by changing the CGHs so that the realized phase pattern closer resembles the desired one [23][24][25][26][27].In this work, however, no such scheme was used.
Besides studying the intensity of the spots directly, we have also observed the effect indirectly by measuring the positions of optically trapped silica spheres, subject to a viscous drag force, as described in section 3.2.The experiment mimics a common application of holographic optical tweezers used in our lab and clearly illuminates the problems that may arise from unwanted intensity fluctuations; an array of traps is to be filled by objects for examination, in our case normally yeast cells, in a microfluidic channel with some fluid flow.Occasionally, an unwanted object is trapped, typically a cell with deviating morphology, which demands for a temporary removal of that spot to release the object.As seen in the previous example, this may cause a decrease in intensity of the remaining spots so that also objects of interest might be lost.
In our experiment, we have trapped a single silica bead in one of the 25 spots arranged in a 5 × 5 grid as previously described and the fluid flow was instead induced by moving the microscope stage.A typical result from such experiment is shown in Fig. 7; in (a)-(c) the speed of the microscope stage was 20 µm/s and in (d)-(e) the speed was 150 µm/s.Even with a relatively small external force, i.e. with low speed of the microscope stage, the displacement of the bead due to the combination of external force and fluctuations in spot intensity is significant in (a), where GS-Rand was used.This may be seen as a sudden deviation from the equilibrium position just after a CGH was sent to the SLM.For increased speed of the  microscope stage, the bead was lost.With both GS-CGH1 and GS-RPC, however, the speed could be increased more than tenfold before the bead was lost.As shown in Fig. 7(d), bead displacement becomes apparent at higher speeds when GS-CGH1 is used while the bead position remains stable during update with GS-RPC.The slow variations in bead position in Fig. 7(d)-7(e) are most likely due to imperfections in the SLM and the optical system that prohibit spot intensities to reach theoretical values as previously discussed.The effect is also clearly seen in Fig. 8(d), 8(h), which shows the intensity of a selection of spots during update between the same CGHs as used in Fig. 7, measured in a free space configuration as described in Section 3.3.Besides lowered uniformity of the intensity in desired spot positions, the intensity in the position of the removed spot is slightly higher than desired.This type of "ghost trap" commonly occurs in symmetric spot arrangements [28] and the effect is more pronounced if the SLM fails to correctly encode the desired CGH.The effect appears randomly and we have not found any general tendency for RPC to increase the problem experimentally.In simulations, however, the ghost traps become slightly more intense for low α due to the restricted design freedom.

Optimizing α
One possible drawback of using GS-RPC may also be observed in Fig. 6; by restricting the maximum change of pixel values, the ability for the algorithm to converge to an optimal CGH may be hampered, resulting in reduced diffraction efficiency.For this reason, the spot intensities after the update has completed are slightly lower in Fig. 6(c) and 6(f) compared to (d) and (e).The severeness of this effect naturally depends on the choice of α; a low value of α means that a larger fraction of the SLM pixels are hindered from reaching a value advantageous for diffraction efficiency.The effect was investigated in a series of simulations where holograms were generated with different values of α for the same sequence of trap positions.The outcome is, however, not only dependent on the choice of α, but on a number of factors such as the number of spots, the step size of the movement and the initial relative phase of the spots.Especially the presence of any symmetry in the spot arrangement was found to affect the ability to suppress intensity fluctuations using GS-RPC.Simulations were therefore performed with different step sizes, with and without symmetry in the spot arrangement and with different number of spots.The good agreement between simulations and experiments justified the use of simulations for this investigation.As generating a hologram and simulating 20 intermediate steps in the transition from the previous one only required approximately 20 ms, a very large number of simulations could be done to obtain statistically significant data.Figure 9 shows the result from the simulations performed in order to investigate the effect of the phase restriction parameter (α) as described in section 4. In each simulation, a number of spots were held stationary while one spot was moved in the x-y plane.Each row of plots corresponds to one type of spot arrangement and movement of spots: in the upper row, spots were positioned and moved randomly; in the middle row, four spots were positioned in a square pattern separated by 10 spot sizes and one spot moved linearly; and in the lower row, 25 spots were positioned on a square grid with 10 spot sizes separation and one spot moved linearly.These three scenarios were intended to cover a broad range of situations common in holographic optical trapping experiments.The use of both random and symmetric arrangements should reveal whether the negative influence symmetry may have on CGHs [28] also influences the behavior during update of the SLM.
The rightmost value in each plot in Fig. 9, where α = 1, correspond to the use of a conventional weighted Gerchberg-Saxton (GS-CGH1) algorithm and is used as a reference in the evaluation of GS-RPC.Not shown in the figure are the results from simulations using GS-Rand.The maximum intensity drop was near 100% in all simulations using GS-Rand and the diffraction efficiency comparable to GS-CGH1.
Each column of plots corresponds to simulations performed with the same step length; the step lengths were (from left to right) 1 spot size, 0.5 spot size and 0.25 spot size.As previously shown in [16], the intensity loss is strongly dependent on the step size and could in principle be reduced to a minimum by using a very small step size.However, this may prove impractical as the limited update rate of nematic SLMs would restrict the maximum movement speed of the spots.Instead, using GS-RPC with a carefully chosen α gives similar reduction of the intensity fluctuations with little or no sacrifice in the quality of the resulting CGHs and maintained maximum speed.For example, using α = 0.88 and steps of 1 spot size gives the same or less intensity loss as reducing the step length by a factor of 4. The corresponding reduction in diffraction efficiency was less than 1% in all examined cases.In general, the previously mentioned reduction in diffraction efficiency with low α was not dramatic; the decrease in diffraction efficiency was less than 1% in all cases where α = 0.90 was used and less than 5% where α = 0.82 was used.The uniformity was found to be arbitrarily close to 1, i.e. perfect uniformity, in all simulations regardless of the choice of α and is therefore not shown in the figure .The third scenario, Fig. 9(g)-9(i), differs slightly from the others in the sense that the simulated maximum increase of the spot intensities is higher instead of lower for high values of α.The reason for this is that the position of one of the spots coincides with the optical axis which is where most of the light is directed during the SLM update.This reflects another issue with the use of traditional algorithms; updating the SLM repeatedly may disturb experiments conducted in the vicinity of the optical axis and even cause damage to light sensitive samples as the temporal increase in light intensity at this position may greatly exceed the intensities of the desired spots.This effect is also reduced by the use of GS-RPC.
We also note that the presence of symmetry in the spot arrangement indeed affects the ability for GS-RPC to suppress intensity fluctuations.While an optimal α may be found around 0.85 if symmetry is avoided [Fig.4(a)-4(c)], a slightly lower α must be chosen to obtain similar improvement with symmetric spot arrangements.The improvement of using α = 0.85 compared to regular GS-CGH1 is nevertheless significant in all cases; we have found 0.85 to be an appropriate choice of α for most situations.

Timing experiments
Many applications of HOT are time critical and much effort has been devoted to minimize the time required for hologram generation [4,5].An important question for the usefulness of GS-RPC is therefore whether our method adds to the computational time of conventional GS.We have performed several timing experiments where identically implemented versions of GS, written in CUDA, with and without the addition of RPC were compared.The time required for hologram generation on our setup was around 7 ms for the FFT based version (nearly independent on the number of spots).For the version based on Fresnel propagation, the computation time was nonlinearly dependent on the number of spots and the total computation time for 3-25 spots is shown in Fig. 10.In the timing experiments, holograms with 512 × 512 pixels were generated in 10 iterations and the presented times refers to complete executions of the program, including transfer of the CGH from the GPU to the host computer.The experiments showed, somewhat surprisingly, that the use of RPC caused virtually no increase in the total computation time.This may be explained by the low arithmetic intensity of the algorithm, i.e. even with the added computations, the execution time is determined by the memory bandwidth of the device, not computational speed.Note that this result is highly hardware specific and only applies to execution on graphic processing units.
As our simple model for estimating the intermediate states of the SLM during update was found to provide accurate results, the simulations may be used to optimize the limiting parameter α as shown in this study.If implemented into the hologram generating program, it could also be used to predict intensity fluctuations in real time during an experiment.This could be useful for example in order to predict unacceptable intensity fluctuations in an OFM experiment.In our implementation, the computation time required for simulation of 10 intermediate holograms was approximately 0.5 ms per spot.

Conclusion
We have demonstrated how unwanted intensity fluctuations occurring during update of an SLM may be virtually eliminated by using a modified hologram generating algorithm.The modification limits the maximum phase difference for each pixel between consecutive holograms and may be implemented into most iterative algorithms.We have also studied how the modification affects the quality of the resulting holograms under various conditions.It was found that in situations where a non-symmetric spot arrangement was used, the decrease in diffraction efficiency was less than 2% while the intensity fluctuations were reduced by more than 95%.The method, which we call RPC, could be implemented into any iterative algorithm for hologram generation.In our experiments, using two different Gerchberg-Saxton based algorithms implemented in CUDA, the modification did not increase the total computational time for the hologram generation.
The technique should be especially useful in optical force measurement experiments using holographic optical tweezers which generally require highly stable and predictable spot intensities.

Fig. 1 .
Fig. 1.Measured far field intensity during update of the SLM (Hamamatsu) with a sequence of ten CGHs (Media 1).Spots were positioned on a 5x5 grid centered on the optical axis.In every second CGH one spot, marked with a circle, at varying positions were removed from the grid and then replaced in the following CGH.Snapshots of the far field are extracted at (a) t = 3.0 s (stable CGH #3), (b) t = 3.26 s (worst case for this specific SLM update), and (c) t = 4.0 s (stable CGH #4), respectively.(d) Examples of normalized spot power as function of time for the specific SLM update shown in (a)-(c).
with an initial guess for the phase distribution, 0 ( , ) x y φ [Fig.2(a)].The complex field in the spot positions is then computed using some numerical propagation method [Fig.2(b)].If the amplitude of the obtained field, ( , , ) n A u v w , does not correspond to some specified requirements, it is replaced by the desired amplitude, ( , , ) n A u v w ′ , while the obtained phase distribution, ( , , ) n u v w φ is maintained [Fig.2(c)].The resulting complex field is then backwards propagated [Fig.2(d)] to the plane of the SLM and the resulting phase

Fig. 2 .
Fig. 2. General flowchart of a GS algorithm utilizing amplitude weighting in the spot positions and a random phase distribution as initial guess for ( , )x y φ

Fig. 3 .
Fig. 3. General flowchart of a GS-RPC algorithm utilizing amplitude weighting in the spot positions and reusing the previous CGH as initial guess for ( , )x y φ

Fig. 4 .
Fig. 4. Schematic of the free space setup described in section 3.1.Lenses L1 and L2 expand the laser beam to match the used area of the SLM (BNS).The far field intensity distribution is imaged onto the CMOS camera using L3.Inset shows a typical CGH used in the experiments.

Fig. 5 .
Fig. 5. Schematic of the HOT setup described in Sections 3.2 and 3.3.Lenses L1 and L2 expand the laser beam to match the width of the SLM (Hamamatsu).The SLM plane is then imaged to the back focal plane of the microscope objective using lenses L3 and L4.The trapping plane coincides with the imaging plane of the microscope.Brightfield images of the trapped bead are captured by a CMOS camera (Prosilica).To allow for the free space experiments described in Section 3.3, a mirror and a second CMOS camera (Photonfocus) was inserted in the optical path as indicated by the dashed box.Inset shows a typical CGH used in the experiments.

y t φ is the phase distribution of an intermediate hologram at time t and 1
distribution of the first and second of two consecutive CGHs, respectively.

Figure 6
Figure 6 shows simulated and experimentally measured intensities of 10 spots during update between two CGHs.This example clearly shows the remarkable improvement of reusing the previous CGH for the calculation of the next (GS-CGH1) compared to the use of a random initial guess (GS-Rand).The maximum intensity drop was reduced by nearly 80% in this case.GS-RPC further reduced fluctuations almost completely.While the simulations are in good agreement with experiments regarding the extreme values of the spot intensities, the

Fig. 6 .
Fig. 6.Spot intensities during SLM (BNS) update of the free space setup described in Section 3.1.The first and second CGHs each forms 9 spots located among 10 equally separated positions on a circle with radius of 70 spot sizes centered around the optical axis; the position of the missing spot being different for the two CGHs.(a) Simulated and (d) measured spot intensities using GS-Rand.(b) Simulated and (e) measured spot intensities using GS-CGH1.(c) Simulated and (f) measured spot intensities using GS-RPC.

Fig. 7 .
Fig. 7. Measured position for a trapped bead exposed to drag force induced by moving the microscope stage at a speed of (a)-(c) 20 µm/s and (d)-(e) 150 µm/s.The dotted lines represent approximate times for the SLM (Hamamatsu) updates.The CGHs were optimized using (a) GS-Rand, (b),(d) GS-CGH1, and (c),(e) GS-RPC.

Fig. 8 .
Fig. 8. Measured far field intensity during update of the SLM (Hamamatsu) with a sequence of ten CGHs optimized using GS-CGH1 (Media 2, (a)-(d) shows data from one update) and GS-RPC (Media 3, (e)-(h) shows data from one update), respectively.Spots were positioned on a 5x5 grid centered on the optical axis.In every second CGH one spot at varying positions, marked with a circle, was removed from the grid and then replaced in the following CGH.Snapshots of the far field are extracted at (a),(e) t = 3.0 s (stable CGH #3), (b),(f) t = 3.26 s (worst case for this specific SLM update), and (c),(g) t = 4.0 s (stable CGH #4), respectively.Normalized spot power for a selection of spots during the specific SLM update shown in (a)-(c) and (e)-(g) are shown in (d) and (h), respectively.Note that the results obtained when using CGHs optimized with GS-Rand is given in Fig. 1.

Fig. 9 .
Fig. 9. Simulated diffraction efficiency after SLM update (black), maximum intensity loss (red) and maximum intensity increase (blue) of stationary spots during update as functions of the restriction parameter (α).In (a)-(c), 4 spots were randomly positioned and one spot moved 1000 steps in random directions.In (d)-(f), 4 spots were positioned on a square grid centered on the optical axis and the lower left spot was moved 100 spot sizes to the left.In (g)-(i), 25 spots were positioned on a square grid centered on the optical axis and the lower left spot was moved 100 spot sizes to the left.The step size was 1 spot size in (a), (d) and (g), 0.5 spot sizes in (b), (e) and (h) and 0.25 spot sizes in (c), (f) and (i).The simulation sequence was repeated 100 times.Each plotted data point is the average along with the standard deviation of all steps and all simulation sequences.

Fig. 10 .
Fig. 10.Total computation time for generation of 512 × 512 pixels holograms in 10 iterations using a Fresnel propagation based version of the Gerchberg-Saxton algorithm implemented in CUDA with and without RPC.