Robust Digital Holography For Ultracold Atom Trapping

We have formulated and experimentally demonstrated an improved algorithm for design of arbitrary two-dimensional holographic traps for ultracold atoms. Our method builds on the best previously available algorithm, MRAF, and improves on it in two ways. First, it allows for creation of holographic atom traps with a well defined background potential. Second, we experimentally show that for creating trapping potentials free of fringing artifacts it is important to go beyond the Fourier approximation in modelling light propagation. To this end, we incorporate full Helmholtz propagation into our calculations.

generic trap shapes, and some shapes have also been successfully realised with laser light 25,26 . However, one practical limitation of MRAF is that in order to achieve good convergence in bright areas of the trapping plane (green areas in Fig. 1), we must surrender control over almost all of the dark, background regions. In particular, the original MRAF calculations provide for only a very narrow region defining the zero background potential to which all the potentials are referenced. This could lead to experimental problems such as inefficient loading of atoms into, or percolation of atoms out of the trap; this issue is particularly important if the desired trap has sharp edges so the region of uncontrolled light intensity is immediately adjacent to the region of high atomic density. If the region of zero potential is extended, the fidelity of MRAF kinoforms is severely compromised. The first result of this paper is to show that this problem can be eliminated by properly redefining the target potential so as to include both the conventionally defined trapping region and a sufficiently large well defined background region (a ''canvas'' on which an arbitrary potential landscape is ''drawn''). We formulate an offset-MRAF (OMRAF) algorithm through a simple set of practical rules which obey the limits imposed by Nyquist's theorem and minimise the occurrence of convergence-stalling optical vortices.
The second result of this paper is to demonstrate that the paraxial approximation often employed in literature 21,27 can lead to undesirable artifacts in the traps. We illustrate the improvement that can be made by direct numerical calculation of the light propagation using the Helmholtz equation. Using our computational algorithms and a simple apparatus, we experimentally demonstrate examples of optical traps suitable for atomic experiments, and compare the results to previous realisations 25,26 .
The paper is organised as follows: First we compare computational results produced using the existing MRAF and our OMRAF algorithms. Then we present and analyse experimental images of light fields shaped by our method. Throughout this investigation, we primarily test our technique by creating an atomtronic OR-gate trap shape 11,21,28 (see Fig 1). To illustrate the universality of our method we also include the final results for a uniform square trap and an annular BEC stirrer 21,25 .

Results
The Gerchberg-Saxton algorithm and MRAF. Within the Gerchberg-Saxton framework, the core algorithm relies on linking the moduli in the SLM and trapping planes by simulating the light propagation back and fourth between these planes. After each propagation, we manually impose the known modulus constraints while leaving the phase to converge on the required solution. The mathematical details are presented in the Methods section, and here we just briefly review the algorithms for kinoform calculation before presenting the numerical and experimental results. It is well documented that convergence of the naïve Gerchberg-Saxton algorithm is highly erratic because the modulus constraints in each plane are nonconvex 29,30 . The MRAF algorithm improves on this by defining a ''drawing'' D in the trapping plane containing all the points in which the desired potential is non-zero, and an additional ''canvas'' region, C, of zero potential around the drawing [ Fig. 2(a)]. The modulus constraint is then manually imposed only inside C and D at each iterative step, thus reducing the number of constraints the algorithm aims to satisfy.
In the original MRAF paper, region C is limited to a tight border around D which has a typical width of around 1% of the trapping region length scales. Although C is a featureless region of zero potential, it is an integral part of the trap since it defines the zero to which all the pixels in D are referenced. For practical purposes such a small canvas can be a severe limitation. Our method aims to expand C to approach its theoretical limit which is set by the Nyquist theorem to be 25% of the area of the simulated trapping plane (see the Methods section). In this large canvas regime, MRAF convergence stagnates and produces poor solutions even after many iterations (see Fig. 2

)
Optical vortices and OMRAF. We observe the problems of the MRAF algorithm for extended C to be due to the formation of large populations of optical vortices during the early iterations (approximately 1 vortex for every 10 pixels in Fig. 2(a)). Optical vortices are topological features of the light field corresponding to a phase winding around a point of zero intensity. Large vortex populations are problematic because individual vortices are difficult to eliminate in further iterations of the algorithm. The reason for this is that the Gerchberg-Saxton algorithm monotonically reduces the RMS error in each iteration 31 , but topological (un)winding operations cause disruption across the entire kinoform. Such global disruption is not in keeping with the monotonic error reduction in later iterations 32 . Therefore any erroneous vortices are ''frozen in'' early on, and the algorithm gets stuck in a local RMS minimum corresponding to a particular vortex distribution.
One possibility for handling vortices is by pairwise creation and annihilation, which requires only local operations; computational algorithms which encourage such pairwise operations are given in Ref. [33]. However, we consider a much simpler solution. Our method simply offsets the trapping pattern inside C and D by a uniform intensity, jDj 2 , to remove all points of zero light intensity. This redefines the zero of the potential but does not change the physics of the trap. With an offset intensity, we know a priori that  we should aim to create no vorticies at all inside C, so we feed the algorithm with an initial kinoform which contains no phase winding. Specifically, we choose a parabolic initial kinoform 21,24 which defocuses the light into a patch with a characteristic size set by C.
The optimal value of the offset intensity is calculated semi-empirically as follows: After just one iteration, the intensity in the trapping plane adopts approximately the correct shape, but exhibits considerable fluctuations around the desired intensity [ Fig. 2(b)]. If any of these fluctuations bring the intensity locally to zero, a vortex may form at this point. Therefore, we increase jDj 2 until the trapping plane contains no vortices inside C after the first iteration [ Fig. 2(b) and (c)], and then trust that very few vortices will be formed in subsequent iterations. For maximum light-usage efficiency, we want to minimise the proportion of the incident light which is steered into the featureless background, i.e. choose the minimal jDj which gives satisfactory results. We find that for all the patterns we tested a suitable offset was jDj , 10 -15% of the maximum amplitude in the trapping plane (i.e. jDj 2 , 1% of the maximum intensity). The computational simulation of the trapping potential produced using the OMRAF algorithm is shown in Fig. 2. Our method reduces the number of vortices seen in C after 30 iterations by a factor of , 20 compared to MRAF. This allows excellent convergence of the algorithm, leaving only 4% RMS error after 30 iterations and 1% RMS after 1000.
One compromise associated with the OMRAF method is a reduction in the efficiency of light use. We define efficiency as the ratio of the integral of the trap intensity above the background offset to the total integral across the trapping plane. Including the non-zero D causes a drop in efficiency by a factor of approximately 2, from 43% to 24%, for the OR-gate. Nevertheless, our method remains more efficient than intensity masking methods such as Ref. [14], which report 3% efficiency for simple patterns. We calculate that for more complex patterns such as the OR-gate, intensity masking would be less than 1% efficient.
Experimental realisation with laser light. We test our OMRAF algorithm experimentally by producing light patterns using 532 nm laser light in the arrangement shown in Fig. 1(a). We use the standard SLM-based Shack-Hartmann algorithm 34 to crudely correct for the low spatial frequency phase aberrations in the optical system and characterise our input laser beam. We then use an active feedback algorithm 25 to optimise the final experimental patterns. This feedback routine is essential to remove the ''sinc envelope'' caused by the pixellation of the SLM, which would otherwise globally modulate the intensity pattern. Finally, in practice, we find that when we try to produce experimental patterns with C covering the maximum theoretically permitted area (25% of the trapping plane), we cannot achieve patterns with less than 20% RMS fluctuations. This is significantly improved in the results presented below by reducing C to cover only 10% of the trapping plane.
In the previous section, we presented computational results generated under the paraxial approximation in which the propagator mapping the light in the SLM plane to the trapping plane is given by a scaled Fourier transform 27 . Using this approach in our experiments with an f 5 200 mm lens, we are able to produce a trapping pattern with RMS error of < 11% [ Fig. 3(a)]. As highlighted in Fig. 3(b), the most prominent form of error is semi-regular fringing. We suggest that the source of this error is the inadequacy of the paraxial approximation, and use a numerical Helmholtz solver (see the Methods section) to test this hypothesis. Specifically, we use the OMRAF algorithm under the paraxial approximation to produce a ''Fourier-generated'' kinoform, and then input this kinoform into our numerical Helmholtz propagator to simulate the experimental apparatus. This indeed produces the fringing pattern similar to the one observed experimentally, which gives us confidence to continue with the Helmholtz method. We thus replace all instances of the Fourier propagator in the OMRAF algorithm with a Helmholtz propagator.
In this way, we are able to eliminate the erroneous fringing, and reduce the RMS error of the experimental pattern to 7% [see Fig. 3(a)]. As shown in Fig. 3(c), we can produce other trap shapes with similar fidelity. Thus, for traps with sharp edges and an extended canvas region, we achieve RMS variation comparable to the 4% fluctuations previously seen only for simple smooth traps 25 . In Ref. [35], it has already been shown that 5% RMS errors are sufficiently low for atomtronic applications. We therefore believe that our methods show great promise for future applications.

Discussion
We have presented a simple method for producing optical traps of arbitrary 2D profile suitable for use in cold atom experiments. A useful optical potential must be referenced to a well defined background, and we showed that in order to do this, we must modify the existing algorithms to include a non-zero background light intensity. Specifically, we developed the OMRAF algorithm guided by the principle that points of zero intensity seed convergence-stalling vortices, and therefore should be avoided. In computer simulations we produce traps covering 25% of the trapping plane (the Nyquist limit) with 4% RMS error after only 30 iterations. We also took the next step of realising these 2D profiles in laser light. Here, we remove systematic aberrations using a Shack-Hartmann method and active feedback. In addition, we found that using a beyond-paraxial Helmholtz solver improves the fidelity of the experimentally produced light patterns.

Methods
Modified Gerchberg-Saxton algorithms. Both the MRAF and the OMRAF approach have the Gerchberg-Saxton algorithm at their core. This algorithm attempts to obtain the kinoform, w(r), which must be applied to the input Gaussian field, E 0 (r), in the SLM plane (spanned by r) to obtain the desired trap shape, T(R), in the trapping plane (spanned by R) at z 5 2f [see Fig. 1(a)]. These two planes are related by a projection operator, P, so we may write the implicit equation for w(r) as: The Gerchberg-Saxton algorithm solves this equation iteratively as illustrated in Fig. 4. MRAF relaxes the modulus constraint outside the trapping region by defining a drawing, D (defined by T(R) . 0) together with a narrow canvas, C, of zero intensity around D, and then applying the following algorithm: We find that a mixing parameter m 5 0.4 gives the minimum RMS deviation for all the patterns trialled (see also Ref. [21]). OMRAF further modifies this algorithm in the following ways: First we expand the canvas so that the trapping region covers the maximum theoretical area set by Nyquist's theorem (25% of the whole trapping plane -see the next section). Secondly, throughout C and D, we offset the desired light intensity according to: This shifts both the trapping potential and its zero reference level equally, and does not change the physics of the trap.
Maximum canvas size. In order to best use the degrees of freedom offered by an N 3 N pixel SLM, we should allow region C to cover as many pixels as possible in the discretized trapping plane. This theoretical limit is found by the following Nyquist argument: The N 3 N SLM plane can be fast-Fourier-transformed to an N 3 N trapping plane with a diffraction limited spot size of exactly 1 pixel. The Nyquist limit states that the diffraction limit should be reduced to half a pixel (the smallest feature in the trap shape). This can be done by padding the computer-simulated SLM plane with zeros to increase its size to 2N 3 2N. This means that only 25% of the simulated SLM plane is actually covered by the physical SLM, so we can at best control only 25% of the trapping plane.
A computationally efficient helmholtz solver. We see improvement in the experimentally realised trapping intensities when the following method is used to model the propagation operator P in the OMRAF iterations: The Helmholtz equation for light propagation in free space gives that the Fourier transform of a light field, E(r; z), is modified by a phase factor H(j; z) as it propagates in the 1z direction 27 where F g r ð Þ; j ½ denotes the Fourier transform (FT) of g, making the Fourier variable pairing r « j explicit, and k 5 2p/l denotes the wave vector of the light. The routine for propagation of E 0 via the apparatus in Fig. 1(a) can then be conceptually decomposed into 3 stages: where represents convolution, and represents the FT of the phase pattern imprinted by the lens, which, after aberration correction, we assume to be: