Adaptive-Additive Algorithm

Our approach is based on the adaptive-additive (AA) algorithm of Soifer et al. [13], an iterative numerical technique which explores the space of degenerate phase profiles, $\Phi^f(\vec \rho)$, to find a phase modulation of the incident laser beam encoding any desired intensity profile in the focal plane. To facilitate calculation and fabrication, both the input and output planes are discretized into $M \times M$ square arrays of pixels. Optimal spatial resolution requires pixels in the focal plane to be one half-wavelength on a side, $\delta^f = \lambda/2$. The number, $M$, of pixels on a side then depends on the desired dimensions of the trapping array. Lengths in the input and focal planes are related by Eqs. (4) and (6), so that the corresponding pixel size in the input plane is $\delta^{in} = \lambda f / (M \delta^f) = 2f/M$. If $\delta^{in}$ is inconveniently small, then L1 and L2 can be chosen so that a more amenable pixel size in the input plane corresponds to $\delta^{in}$ in the focal plane.

Figure 3: Flow chart for the adaptive-additive algorithm. The phase modulation, $\Phi^{in}_n(\vec r)$, can be quantized into discrete steps with every iteration, as shown, or after the algorithm has converged.

The AA algorithm, depicted in Fig. 3, starts with an arbitrary initial guess for $\Phi^{in}_1(\vec r)$ and an initial input wavefront $E^{in}_1(\vec r) = E_0(\vec r) \, \exp[ i \Phi^{in}_1(\vec r)]$. The Fourier transform of this wavefront is the starting estimate for the output electric field: $E^f_1(\vec \rho) = {\cal F}\{E^{in}_1(\vec r)\} = A^f_1(\vec \rho) \, \exp[i\Phi^f_1(\vec \rho)]$. The corresponding intensity in the output plane, $I^f_1(\vec \rho) = \vert A^f_1(\vec \rho)\vert^2$ is unlikely to be a good rendition of the desired intensity pattern, $I^f(\vec \rho) = \vert A^f(\vec \rho)\vert^2$. The error,

$$\epsilon_1 \equiv \frac{1}{M^2} \sum_{i=1}^{M^2} [I^f(\vec \rho_i)-I^f_1(\vec \rho_i)]^2,$$ (16)

is reduced by mixing a proportion, $a$, of the desired amplitude into the field in the focal plane:

$$\bar E^f_1(\vec \rho) = [ a A^f(\vec \rho) + (1-a)A^f_1(\vec \rho) ] \, \exp [i \Phi^f_1(\vec \rho)].$$ (17)

Inverse transforming $\bar E^f_1(\vec \rho)$ yields the corresponding field in the input plane, $\bar E^{in}_1(\vec r) = \bar A^{in}_1(\vec r) \, \exp[i \Phi^{in}_2(\vec r)]$. At this point, the amplitude in the input plane no longer matches the actual laser profile, so we replace $\bar A^{in}_1(\vec r)$ with $A_0(\vec r)$. The result is an improved estimate for the input field: $E^{in}_2(\vec r) = E_0(\vec r) \,\exp[i \Phi^{in}_2(\vec r)]$. This completes one iteration of the AA algorithm. Subsequent iterations lead to monotonically improving estimates, $\Phi^{in}_n(\vec r)$, for the desired phase modulation [13]. The cycle is repeated until the error, $\epsilon_n$, in the $n$-th iteration converges to within an acceptable tolerance: $(\epsilon_n - \epsilon_{n-1}) / \epsilon_n < \chi$.

The phase and amplitude fields are computed as arrays of double-precision numbers, and their Fourier transforms calculated with fast Fourier transform (FFT) routines. Starting from random input phases, $\Phi^{in}_1(\vec r_i)$, uniformly distributed in the range 0 to $2 \pi$, the AA algorithm typically requires eight iterations to converge within $\chi = 10^{-6}$ of an acceptably accurate local minimum of $\epsilon_n$ using an intermediate value for the mixing parameter, $a = 0.5$.