1. Introduction

Computer generated holograms (CGHs) for phase modulation have become an established optical technology for efficiently producing an almost arbitrary spatial intensity distribution using a single optical component. In particular CGHs have proven very powerful in generating focused spots in an arbitrary 3D configuration in space – a function which can be regarded as sophisticated beam splitting, with arbitrary directions, intensities and divergence (or convergence) of the individual beamlets. This is the type of CGH considered in this work. Traditionally, these CGHs, also known as diffractive optical elements or kinoforms, were tailor-made with costly lithographic techniques[1, 2]. With the rapid development of highly accurate phase-only spatial light modulators (SLMs) it has become possible to realize sequences of CGHs dynamically. The physical phase modulation is then obtained after the response time of the SLM, which is typically of the order of milliseconds. Thus, spot-generating CGHs produced by SLMs are increasingly being utilizied in applications such as optical interconnects[3], beam steering[4–6], and optical tweezers[7–10]. With ever faster SLMs, the main obstacle in real-time applications is the time needed to calculate the phase-modulating function that solves the desired design problem, i.e., that yields an CGH which after being displayed on the SLM produces focal intensities of desired strength at the desired locations – and none elsewhere. For applications is optical tweezers it would ideally be possible to calculate a phase function for new design task at rates of ~20–50 Hz. For state-of-the-art SLMs containing several hundred thousand individually controlled pixels, this is an immense task, since each pixel needs to be accounted for in the numerical design method. The method we have developed is an effort to reduce the computational time needed when calculating highly efficient CGHs capable of performing complex beam positioning tasks.

For an accurate, efficient, and controlled spot generation, SLM-optimizing design methods should be utilized. Two types have traditionally been widely used. The first type is optimization based on direct search[11–13], in which typically all pixels are sequentially tested for a set of possible values for the phase modulation. For each pixel, a new value is selected based on its ability to increase (or decrease) some quality/merit (cost) function. After all pixels have been optimized once, the procedure normally needs to be repeated a few times since the optimal phase modulation of one pixel is influenced by the phase modulation of all other pixels. To avoid stagnation in a local optimum, some degree of randomness can be included in selecting the optimal phase value for each pixel, for instance with simulated annealing[14–16]. Although always feasible in theory, direct search methods are very computationally demanding, and thus impractical in near real-time applications.

The second type of design method is based on the Gerchberg-Saxton (GS) algorithm[17]. Rather than optimizing the phase of individual pixels, this method is based on projecting an entire field onto a field with the desired properties in a single step. In practice, the projections in the GS algorithm are simple changes of numerical values in the calculated optical fields. At the SLM the calculated field is changed so that it can be generated by the SLM, e.g., so that it contains only phase modulation. The calculated diffracted field, on the other hand, is changed to better agree with the desired, according to some specifications. These projections are done sequentially, alternating between the SLM field and the diffracted field with a propagation calculation in between. Each projection partly destroys the effect of the proceeding projection and thus this is an iterative process that is repeated until the resulting fields are no longer changed significantly by the projections. However, since in the projection method there is no need to account for the contribution of an individual pixel to the diffracted field, as in a pixel-by-pixel optimization, the GS-based design is considerably faster than the direct search methods. There exist a large number of modified GS algorithms, or similar projection-based optimization methods, for the design of spot generating CGHs[18–22].

Having chosen the GS-based optimization because of its huge speed advantage, there is yet another decision to be made that influences the speed – the choice of the numerical propagation method for the back and forward propagation between the SLM and the evaluation plane(s) for the diffracted field. Each iteration involves at least one forward and one backward propagation calculation, so having an efficient propagation algorithm is crucial for maximum speed. Most SLMs have pixel sizes much larger than the wavelength of the light and thus the SLM produces a highly paraxial field whose propagation can be excellently treated within the Fresnel approximation[23]. This simplifies the mathematical treatment to essentially a Fourier transform, or its inverse for the back-propagation.

The numerical Fourier transform can be calculated straightforwardly such that the field is obtained only in those positions where the focal spots should be positioned, and nowhere else. In this case, the spot positions can be arranged in 3D, and they can also be individually shaped, e.g., as ordinary Gaussian or Laguerre-Gaussian beams of different orders[19]. The straightforward calculation of the Fourier transform makes the method fairly slow, particularly for larger number of focus positions, since the computation time increases linearly with this number.

A different type of Fourier transform is the fast Fourier transform (FFT), which is the most commonly used in GS-based methods. Due to a highly clever calculation method, the Fourier transform is calculated in a very large number of positions – N×N, which is the same as the number of sampling positions for the entire field in the plane of the SLM – in a very short time, although not quite as short as the calculation time for a single position using the direct method. However, the high calculation speed comes at a price; even though the steering accuracy of SLMs inherently is very high[25], one must accept that with the FFT the Fourier transform is calculated on a rectangular equidistant grid of positions located in a single plane perpendicular to the direction of propagation. In addition to the speed advantage, the large number of evaluation positions for the FFT implies that in principle the entire field emitted from the SLM can be monitored and manipulated at a certain distance. Thus, unlike the direct method, at only little extra computational cost one can avoid unintentional intensities at positions where there should be no focal points. Evidently, such a functionality is crucial for the proper operation of optical tweezers. In spite of this, the FFT might seem disappointingly wasteful, since generally very few of the evaluation points for the Fourier transform are actually used. It may also seem overly restricted, due to the requirements on the location of the evaluation points. Indeed, it is possible to obtain a 3D focal arrangement, but this is done by introducing a new calculation plane, and one additional FFT, for each unique focal distance from the SLM, which drastically increases the calculation time[26–28]. Moreover, even though it is possible to form different images in the different planes, it is not possible to control field characteristics smaller than the grid point separation, e.g., to let some of the spots have a Laguerre-Gaussian shape.

In this work, we have developed a single-plane FFT-based design method that largely circumvents these drawbacks. By using a single plane, the iteration plane (IP), which does not in general coincide with the location of any of the focal positions, grid-free positioning in 2D or 3D is possible with an approach similar to the conventional GS-based design. The only restriction is that in the IP there should be only a small spatial overlap between light going to different focal positions. This means that the method is well suited for 3D problems with focal positions that are clearly separated transverse to the direction of light propagation, and whose separation in the longitudinal direction, along the light propagation, is not too large. The grid-free 3D positioning is made possible by using more of the information contained in the single IP, resulting in a design method that is almost as fast as the simplest and most restricted – and most rapid – implementation of the GS algorithm. In fact, the new design method deviates significantly from the simplest single-plane GS method in only two respects: a slightly different projection in the IP, and a slightly more elaborate procedure to find the desired field in this plane (toward which the propagated field is then being projected). This increases the time required for the initiation part, before the iterative part of the design, but it can be reduced to constitute only a marginal fraction of the total execution time, as will be explained in the following. We will also explain the other relevant details of our method, and show simulated and experimental results that agree very well with the different design tasks.

 

Fig. 1. (a) An SLM and its output field consisting of two spots positioned at different distances from the SLM. The iteration plane (IP) and the optimization regions are indicated. Typical (b) amplitude and (c) phase of E^desired_n for an optimization region.

Download Full Size | PDF

2. The FFT-based design algorithm

Our novel FFT-based algorithm is based on the simple fact that the optical field distribution giving rise to a focused spot with a specified intensity, phase, and location, is known at any distance from the spot itself. In Fig. 1(a), the optical field propagating from the SLM is shaped such that two spots are formed at two different distances from the SLM. The basic principle of our approach is thus that in the IP, the optical field in each of the two indicated regions is given by the position and the complex amplitude of the corresponding spot. An example of the amplitude and phase of an optical field in the IP, which produces a spot at some distance from the IP, is shown in Figs. 1(b) and 1(c), respectively. The only thing needed for the design algorithm to determine is a phase offset, i.e., the addition of a constant phase value to the field in each optimization region. This does not change the spot intensity or position. The phase offset values for the optimization regions in the IP are thus the optimization parameters in our algorithm. This is in exact correspondence with the conventional GS-algorithm, which uses the phase of the field in each spot position as the optimization parameter. Evidently, an addition of a constant phase to the field in an optimization region in the IP changes the phase in the corresponding spot position by the same amount.

The flowchart of the algorithm is shown in Fig. 2(a). The initializing step of the algorithm calculates the desired optical field (amplitude and phase distribution) in each optimization region n of the IP. We denote the amplitude distribution by A ^desired _n and the phase distribution by φ ^{desired, no offset} _n since we define this initially calculated phase distribution to have the phase offset value ϕ_n=0. The iterative part of the algorithm starts by propagating the so obtained desired field in the IP, E ^desired _IP, back to the SLM plane, resulting in E ^desired _SLM. In the SLM plane, phase freedom applies just as for the conventional GS-algorithm, i.e., the SLM is assumed to be able to realize any phase value, and thus the desired phase distribution is left unchanged. However, the phase-only SLM cannot affect the amplitude distribution which is hence set to that of the incident laser beam.

After propagating this field to the IP, resulting in the field E ^obtained _IP, the power P ^obtained _n in each of the optimization regions n is calculated as well as the diffraction efficiency $\eta =\sum P_n^{\mathrm{obtained}}/P_{\mathrm{total}}$ and the uniformity error $\epsilon =\genfrac{}{}{0.1ex}{}{\max (p_n)-\min (p_n)}{\max (p_n)+\min (p_n)};p_n=P_n^{\mathrm{obtained}}/P_n^{\mathrm{desired}}$ where P _total is the total power in the IP and P ^desired _n is the desired power in optimization region n.

 

Fig. 2. (a) Flowchart of the algorithm. (b) Phase read-out and reshaping procedure for a single optimization region. (c) Amplitude modification scheme for the IP-field. Note that the amplitude is shown in false scale to enhance the visibility of the unwanted ghost spots.

Download Full Size | PDF

In the IP, the phase offset value ϕ_n is determined for each optimization region. The phase offset value should represent the “phase difference” between the obtained field in the optimization region and the field from the initializing step. This entity can be extracted, e.g., as a weighted average of the phase differences in the pixels included in the optimization region n, with A ^desired _n as the weighting function. However, we found that it was sufficient to use the phase difference in the “most significant pixel”, i.e., the one which has the highest value of A ^desired _n. After having found the ϕ_n-values from this iteration, they are added to the initial phase distributions as indicated in Fig. 2(b).

As mentioned, the desired amplitude distribution in the IP is already determined by the amplitude distribution in the optimization regions, A ^desired _n, from the initialization step. Elsewhere in the IP it should be zero. However, as in many GS-based design methods the projection optimization can be made more efficient by adding some freedom for the desired amplitude in each iteration. Therefore, to improve the uniformity among the desired spots, a weighting procedure is used, by which the desired amplitude in optimization region n is multiplied by a weight w_n. The weights work to boost and supress the optical power in regions that tend to contain too little and too much power, respectively. The weighting process can be done in many ways[18, 19, 29–31]. Here, we use the weighting defined in [31] as

The convergence parameter c, which generally is chosen to take values between 0.05 and 0.3, was set to 0.25 for the optimizations shown in this work.

In order to further increase the uniformity among the different spots and to speed-up the convergence of the algorithm we also use a peripheral region of no interest (RONI) in the IP in which some amplitude freedom applies[32, 33]. In our case, a circular region of interest (ROI) which is centered on the optical axis and covers all desired optimization regions is used, see Fig. 2(c). The area outside the ROI forms the RONI. As long as the power contained in the RONI does not exceed some small predetermined fraction of the total power in the IP, the field in the RONI is not suppressed in the iterations. Should the power be too large, the field is only suppressed to the point where the power is the maximum allowed. The drawback with this amplitude freedom is evidently that some additional power is lost to the RONI, but on the other hand it is possible to direct undesired ghost spots to fall into the RONI, yielding a higher signal-to-noise (SNR) ratio within the ROI. The complete procedure of modifying the amplitude of the obtained IP-field is illustrated in Fig. 2(c). The finally obtained IP is then propagated back to the SLM plane again, and the iterations continue until the diffraction efficiency and uniformity error in the IP, generated by the SLM, fulfil the design criteria.

2.1. Determination of the desired E-fields

A cornerstone of the described method is the initialization step, in which the position and size in the IP of each optimization region, and its desired field E ^desired _n, are determined. Roughly, the size and position of each optimization region can be obtained from the geometrical construction shown in Fig. 1(a) once the position of the IP is established. This position is not very critical but should not be too far away from the focal positions of the spots, since the optimization regions should be small enough to avoid overlap. In fact, it is sufficent to use optimization regions of only a few pixels in side length to ensure grid-free high-precision 3D positioning of the corresponding spots. Thus, for a 2D (or near-2D) desired configuration of spots, the IP may be positioned shortly before or after the spot positions, while for a more extended 3D configuration it is advantageous to position the IP near the center of the region where the spots should be located, such as the case shown in Fig. 1(a).

The most direct method to determine the desired field in an optimization region, E ^desired _n, is to create the spherical wavefront in the SLM plane that has its center of curvature at the desired spot position and then propagate this field to the IP by the same FFT-based propagation method that is used in the iterative part of the algorithm. To avoid calculating one FFT operation for each desired spot, the desired SLM field is calculated as the complex sum of all the spherical wavefronts corresponding to all the desired spots. This complex field (note that no amplitude restriction on the field in the SLM plane is used in this initialization step) is then propagated to the IP with a single FFT operation. For each optimization region, E ^desired _n is then obtained as the portion of the field in the entire IP that is contained in the particular optimization region.

 

Fig. 3. Schematic free-space setup. Lenses L₁ and L₂ expand the HeNe laser beam. Polarizer P is oriented such that the SLM yields phase-only modulation. The intensity distribution at some distance from the SLM is captured with a CCD camera. Inset shows the 256×256 pixels of a calculated CGH (phase pattern) displayed on the SLM.

Download Full Size | PDF

An even more rapid approach, which we have used exclusively in this work, is to use the analytical propagation formulas for Gaussian beams. Then, E ^desired _n is directly calculated from the lateral displacement of the corresponding spot from the optical axis and the distance between the spot and the IP. To allow for creation of Laguerre-Gaussian shaped spots, a generalized expression for Laguerre-Gaussian modes was used[24]. The advantage of this method is that it is less computationally demanding and thus is very fast. The drawback is that the method works best in cases when the incident laser beam has an approximately Gaussian shaped intensity distribution.

After the desired field, E ^desired _n, has been calculated, the actual size of the optimization region is finally determined such that only pixels in which the value of A ^desired _n is a significant fraction of the maximum value of A ^desired _n are included in optimization region n. This avoids the use of unnecessarily large optimization regions and thus minimizes the risk of spatial overlap between neighboring optimization regions.

3. Performance of the algorithm and the obtained SLM frames

To verify the capabilities of the algorithm, CGHs were optimized for different desired output intensity distributions. The evolution of the algorithm was evaluated by studying the diffraction efficiency η and the uniformity error ε with the number of iterations.

We used N=256 and wavelength λ=543 nm. The pixel separation of the SLM was d _SLM=15 µm. A maximum of 10% of the total optical power was allowed to fall outside the ROI. Even though good results typically are obtained after ~15–25 iterations, 50 iterations were used in all our designs.

The simulated intensity distributions from the obtained CGHs were analyzed in detail. To verify that the calculated phase patterns could be realized as SLM frames, measurements using a commercial SLM (Boulder Nonlinear Systems 512×512 SLM system) were performed. A schematic of the free-space setup is shown in Fig. 3. As a light source, a HeNe laser (λ=543 nm) was used. To make sure that the SLM frames were encoded correctly, a diffraction based characterization method[34] was used to measure the phase modulation characteristics of the SLM. At the used wavelength, the SLM allowed for phase-only modulation over a range of [0,2π] radians. To reduce undesired effects caused by the peripheral part of the slightly curved SLM backplane, only the central 256×256 pixels were used in the experiments. The other pixels were set to form a binary grating deflecting the light incident on these pixels to positions outside the region of interest for the measurements.

 

Fig. 4. (a) The desired 5×5 spot pattern with movement directions of the four moving spots indicated. (b) The evolution of the efficiency η and the uniformity error ε of the power in the IP and SP for the design of one of the 25 CGHs. Note that the uniformity error in the IP and SP are close to identical. Simulated (c) and measured (d) relative positions for the 25 spots of the 25 CGHs.

Download Full Size | PDF

3.1. Grid-free positioning in 2D

To verify that the positioning of spots is not limited to a grid, 25 CGHs each producing 5×5 spots in a plane at a distance z=0.50 m from the SLM were optimized. The IP was positioned at a distance zIP=0.46 m from the SLM. For the first CGH the central spot coincided with the optic axis and the spots were regularly spaced with a separation of 1.17 mm, see Fig. 4(a). The coordinates of four of the 25 spots were changed in the CGH sequence, moving the four spots in nominal steps as small as 1/25 of the grid point spacing d _FFT=λz/(Nd _SLM). A typical evolution of the CGH performance during the design with the algorithm is shown in Fig. 4(b). Note that the agreement between the quantities calculated in the IP and in the actual spot plane (the latter calculated only for monitoring the algorithm performance) is very good.

The designed CGHs were evaluated numerically by calculating highly spatially resolved fields around the spots and obtaining precise spot locations by center-of-mass calculations. In Fig. 4(c) these results are shown as the deflections from the start positions. The CGHs were then realized on the SLM and the spot positions were determined in a similar way with center-of-mass calculations based on the intensity distribution captured by the CCD camera. The measured results are shown in Fig. 4(d). The simulated and measured radial position accuracy were within ±0.10dFFT and ±0.12d _FFT, respectively. Calculations and measurements on another series of CGHs showed that the simulated and measured radial position accuracy was ±0.033d_FFT and ±0.035d_FFT, respectively, if only one of the 25 spots was moved. Thus, the agreement between simulations and measurements is very good.

 

Fig. 5. Desired (a) amplitude and (b) phase in the IP (positioned 46 cm from the SLM) for Laguerre-Gaussian of (0,1) order (outer circle) and Gaussian (inner circle) shaped spots, focused in a plane 50 cm from the SLM. (c) The evolution of the efficiency η and the uniformity error ε calculated from the power distribution in the IP and SP. (d) Spatial phase modulation of the optimized CGH. (e) Simulated and (f) measured intensity distribution in the spot plane.

Download Full Size | PDF

3.2. Laguerre-Gaussian shaped spots

We also verified that the algorithm allows different spots to take on different shapes. The desired output intensity pattern consisted of two concentric circular arrangements of differently shaped spots, see Fig. 5. The inner and outer circle consisted of 9 Gaussian and 12 Laguerre-Gaussian (LG) of (0,1) order shaped spots, respectively. The plane of the spots and the IP were positioned z=50 cm and zIP=46 cm from the SLM, respectively. The optimization regions had sizes of 9×9 and 11×11 pixels for the Gaussian and LG shaped spots, respectively.

As seen in Fig. 5(c), the agreement between the uniformity error in the IP and the SP is very good. However, the efficiency figure is ~8% lower in the SP than in the IP since we could use narrower boundaries for the extension of each spot in the SP, thus making sure that only the light actually going to the spots was accounted for as useful in the efficiency calculations. The simulated and measured output field clearly resembles the desired pattern in the SP, see Figs. 5(e) and 5(f). The deviation from circular symmetry of the LG spots in Fig. 5(f) is caused not by the algorithm, but by the slightly curved back plane of the used SLM, which induces a slight astigmatism.

 

Fig. 6. (a) Schematic 3D design task (not to scale); 25 spots positioned in groups of five at z ₁=45.0 to z ₅=55.0 cm using z _IP=49.0 cm. (b) Amplitude and (c) phase of the desired field in the IP. (d) The evolution of the efficiency η and the uniformity error ε based on the power distribution in the IP and the power in the actual spot positions. (e) Measured intensity distributions in the planes, parallel to the SLM, containing the spots. Note that the intensity is shown in a logarithmic scale to enhance the visibility of the defocused spots.

Download Full Size | PDF

3.3. Spot positioning in 3D

To verify that the algorithm allows for 3D positioning of spots, it was used to solve the design task schematically shown in Fig. 6(a). As seen, 25 spots were designed to be grouped in five lines positioned in five different planes at z ₁=45.0, z ₂=47.5, z ₃=50.0, z ₄=52.5, and z ₅=55.0 cm from the SLM. If geometrically projected onto a plane, the focused spots would form a regular square pattern with a spot separation of 1.17 mm. The IP was positioned at z _IP=49.0 cm. The square shaped optimization regions had a side length of 9, 3, 3, 9, and 13 pixels for the plane positions z ₁ to z ₅, respectively, see Figs. 6(b) and 6(c).

As can be seen in Fig. 6(d), the convergence behavior is similar to that of the other design tasks, again with a good agreement between the entities calculated in the IP and the SP. Finally, the resulting CGH was displayed on the SLM and the resulting intensity distribution was captured by the CCD camera in the five planes containing spots. As seen in Fig. 6(e), which shows the captured intensity distributions in a logarithmic scale, the algorithm clearly allows for 3D positioning of spots using a single iteration plane.

4. Conclusions

We have presented an FFT-based design method which optimizes a single-plane defocused field corresponding to a desired spatial configuration of focused spots. Since much more of the information-carrying capacity of the sampled propagated field is utilized in the optimization, design tasks which are difficult with conventional FFT-based design methods can be solved at virtually no extra cost compared to conventional design for standard design tasks. The proposed design method is based on a projection optimization, similar to that in the commonly used Gerchberg-Saxton algorithm. The main difference is a new initialization step and a slightly more complex update operation in the iteration plane.

The capabilities of the new algorithm were verified with three design tasks particularly difficult for conventional FFT-based design methods. First, the algorithm was shown to allow for grid-free positioning in 2D. The simulated and measured smallest position change for a spot was only one eighth (or smaller) of the smallest possible position adjustment with the spots confined to the FFT grid. Second, we verified that a phase hologram producing a mixture of Laguerre-Gaussian shaped spots could be designed. Finally, the algorithm was used for 3D design of a phase hologram which positioned 25 spots, in groups of five, in five different planes. Also for the two latter cases, the measured intensity distributions to a high accuracy agreed with the desired field distributions.