Experimental demonstration of holographic three-dimensional light shaping using a Gerchberg–Saxton algorithm

We use a three-dimensional Gerchberg–Saxton algorithm (Shabtay (2003) Opt. Commun. 226 33) to calculate the Fourier-space representation of physically realizable light beams with arbitrarily shaped three-dimensional intensity distributions. From this representation we extract a phase-hologram pattern that allows us to create such light beams experimentally. We show several examples of experimentally shaped light beams.


Introduction
A monochromatic light beam propagating in free space is a three-dimensional (3D) field and its intensity distribution forms 3D patterns. At the same time, the field cross-section in any one transverse plane completely determines the field everywhere else. This dependence of the entire beam on its cross-section in a single plane forms the basis of beam propagation methods (see for example [1]). It also restricts the possible 3D shapes a monochromatic light beam's intensity can take.
Many methods allow shaping of light beams in 2D, ranging from passing the beam through a shaped aperture to algorithms for calculating phase-hologram patterns that shape the beam in a subsequent plane (see for example [2]). Shaping light beams into non-trivial but restricted 3D shapes is possible by utilizing phenomena such as the Talbot effect [3] or spiral-type beams [4]. It is also possible to shape a light beam into 3D configurations of bright spots [5], whereby the spots can be shaped individually [6]; these techniques are important in the field of holographic optical tweezers. The 3D beam shaping into arbitrary shapes has been demonstrated using computationally intensive direct search methods 2 [7], but continues to be a challenge.
A recent algorithm [8] finds beams whose shape is an approximation within physically realizable limits to any arbitrary 3D target shape. We apply this algorithm to calculate examples of light beams shaped in 3D. From the beams' Fourier-space representations we calculate phaseonly hologram patterns that allow us to create the beams experimentally.
This paper is organized as follows. Section 2 reviews the Gerchberg-Saxton (GS) algorithm [2], which allows light shaping in one plane and forms the basis of what is to follow. Section 3 discusses the 3D generalization of the GS algorithm [8] that we use to calculate shaped light beams. Some of the details of our implementation of the 3DGS algorithm can be found in section 4. Finally, section 5 describes our experiment for the generation of shaped light.

GS algorithm
The GS algorithm [2] is an iterative method originally developed for recovering the phase of an electron or light beam from its intensity distributions in two transverse planes. It can also be applied to shape a light beam, specifically to calculate the phase pattern which light at one plane would require to form an almost diffraction-limited approximation to any desired intensity pattern at a second plane. As the phase pattern required in the first plane can be imprinted onto the light beam with a phase hologram, this allows 2D holographic light shaping. If the requirements in the two planes cannot be met simultaneously within what is allowed by the laws of diffraction, the GS algorithm finds a useful compromise to reconcile these conflicting requirements [9]. For numerical simplicity, the second plane is usually chosen to be the far field of the first plane; mathematically, the fields in the two planes are then Fourier transforms of each other. Here we restrict ourselves to this case.
We write the intensity distribution of the unshaped light beam in the first plane-the hologram plane-as I H and the target intensity distribution in the Fourier plane as I T . Both I H and I T are functions of x and y, which are stored in the computer as 2D arrays of real-valued numbers. The light fields in the two planes are represented by 2D arrays of complex numbers, 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT u n in the hologram plane andũ n in the Fourier plane; the subscript n indicates the iteration number. The initial phase distribution in the hologram plane, ϕ H 0 , is set to any arbitrary distribution, often uniform, i.e.
One iteration of the algorithm, which calculates an improved phase distribution in the hologram plane, ϕ H n , from the previously calculated phase distribution, ϕ H n−1 , then progresses as follows: (3) Equation (2) calculates the initial field in the hologram plane, from which equation (3) then calculates the phase distribution in the target plane, ϕ T n . It uses a fast Fourier transform (FFT), or more generally any discrete Fourier transform. Equation (4) combines this phase distribution in the target plane with the target intensity, I T , giving the field u T n , from which equation (5) then calculates the corresponding phase distribution in the hologram plane, ϕ H n . Over a number of iterations, the actual intensity in the target plane, I n = |u T n | 2 , converges to an almost diffractionlimited approximation of the desired intensity there, I T ; ϕ H n is the corresponding phase-hologram pattern needed to produce this field in the target plane.

GS algorithm in three dimensions
The fields in the two planes in the GS algorithm can also be seen as a field in a real-space plane (x, y) and its k-space (k x , k y ) representation. This allowed Shabtay [8] to adapt the GS algorithm to 3D: instead of dealing with fields in two planes, the algorithm deals with fields in two volumes. One field is a 3D real-space representation of the part of the beam that is to be shaped, the other is its 3D k-space representation, which needs to be consistent with the beam's wave number spectrum. The two fields are 3D Fourier-transform pairs (figure 1). Note that the 2D GS algorithm usually shapes light in the Fourier plane, while the real-space plane contains restrictions such as the intensity profile of the unshaped beam, whereas the 3D GS algorithm shapes light in the real-space volume while the Fourier volume contains the physical constraints.
As discussed in section 1, the field in a single plane completely describes a monochromatic light beam. This is reflected by the fact that the k-space representation of such a beam is confined to a surface. In monochromatic light of wavelength λ, the wave-vector component in the direction of propagation, k z , is related to the values of k x and k y through the equation where k 0 = 2π/λ. This forms a sphere of radius k 0 in 3D k space-an Ewald sphere [8]. We consider here light moving in the positive z direction only, hence the wave vectors are limited to the half of the sphere with k z > 0 (figure 1).  We have written a C ++ implementation of the 3D GS algorithm. Figure 2 shows examples of light shaped using this program. The algorithm can clearly produce intensity distributions in which the target shapes can be recognized. When, as in figure 2, the intensity is shown inverted, the intensity distributions sometimes look like some form of 3D expressionist charcoal drawings.

5
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Details of our implementation of the 3D GS algorithm
In this section we explore some of the details of the implementation of the algorithm. We use the error [10] x,y,z I T (7) to quantify the similarity between the target intensity distribution, I T (x, y, z), and the intensity distribution the algorithm produces after n iterations, I n (x, y, z).

Thickness of the k-space sphere
A purely monochromatic beam restricts the allowed k-space elements to a spherical surface. In our numerical simulations, the 3D k space is discretized, and it is not immediately obvious how to represent the spherical surface that corresponds to a monochromatic light beam in this discrete space. At the same time, the beam's real-space representation is restricted to a cuboid volume of finite size, which implies that the light's wavelength cannot be represented exactly.
In an effort to understand this better we investigate here one aspect associated with the discrete nature of a light beam whose 3D light field and k-space distribution are represented on cubic grids of points. Figure 3 shows the radial cross-section of the discrete-k-space representation of a simulated monochromatic light beam. This k-space representation was calculated by taking the 3D discrete Fourier transform of the beam's spatial representation over a discrete 3D volume. The spatial representation was calculated using a standard beampropagation algorithm 3 [1], starting with the field cross-section of a Gaussian beam close to the beam waist and propagating it into 63 further, equally-spaced, planes. In addition to a peak centred at k 0 = 2π/λ, the radial profile of the resulting discrete-k-space-hemisphere distribution has other distinct features. We believe these additional features to be due to the beam being represented only in a cube and the represented beam therefore having a top-hat profile in the x, y and z directions, which in turn leads to the k-space distribution being widened, more precisely being convolved with a sinc function in the k x , k y and k z directions.
In our simulations, we use perhaps the simplest form of the k-space hemisphere: for each represented pair of k x and k y values, the power in the discrete k z value closest to (k 2 0 − k 2 x − k 2 y ) 1/2 is set to one and that of all the other k z values to zero. Such a k-space hemisphere is exactly one element thick in the k z direction; the radial power profile would be a top hat of width 1 element. This is clearly different from the arguably 'natural' representation shown in figure 3, but it works well enough to produce good experimental results (see below) and simplifies the extraction of phase-hologram patterns (see section 5.1). However, it is clearly an aspect of our implementation that could benefit from further investigation.
It is worth discussing very briefly the case of thicker hemispheres. Figure 4 shows results from the 3D GS algorithm with a k-space hemisphere with a Gaussian radial power profile of variable width. It can be seen that an increase in the thickness of the hemisphere leads The plot shows the powers in all discrete-k-space elements with 0.9k 0 k r 1.1k 0 (where k 0 = 2π/λ) as a function of the radial wave number, k r , calculated for a tightly focused monochromatic Gaussian light beam. The discrete k-space representation was found by numerically calculating the light beam on a 64 × 64 × 64 grid representing a cube of side length 25λ (where λ is the wavelength of light), which was then 3D-Fourier transformed. Whereas in continuous k space all the power would be at k r = k 0 , in discrete k space the power is distributed within a few k-space-element widths of k 0 . Because the edges of the represented cube act like hard-edged slit apertures in the x, y and z directions, the light field's k-space representation is that of continuous k-space, convolved with sinc functions in the k x , k y and k z directions (it is, of course, also discrete).
As the beam used in this example propagates mainly in the positive z direction, the radial k-space structure of the power density broadly reflects the sinc 2 -type structure in the k z direction.
to lower errors. The reason for this is that a thicker k-space hemisphere implies more nonzero k-space elements whose phase the algorithm can alter. Physically, a k-space hemisphere of non-zero thickness corresponds to the presence of a range of wavelengths 4 and therefore polychromatic light; the represented field is a snapshot of the light field at one particular instant. At this instant, the plane-wave components of all colours have the relative phase given in the k-space representation; as different colours correspond to different frequencies, different colours where k 0 = 2π/λ is the radius of the hemisphere; the power of k-space elements with k z 0 is set to zero. The light beams corresponding to some of the points are also shown, each from two different perspectives; the propagation direction (red arrow) is always in the vertical direction. Also shown are the target shape (top inset on the right) and the light beam resulting from a single-element top-hat radial profile (lower inset on the right), which has an error comparable to a τ = 0.0041, which for the grid size used (128 × 128 × 128) corresponds to a Gaussian full width of approximately 1 element.
change their relative phase immediately afterwards. The time evolution of the instantaneous field in polychromatic light can be controlled by an extended GS algorithm that incorporates a time dimension added to the real-space representation of the field and a frequency dimension added to the k-space representation. Experimental realization requires that the relative phase between the plane-wave components of all colours can be controlled. This can be seen as an extension of the shaping of the time-resolved field of short pulses, which has previously been demonstrated experimentally [12].

Influence of numerical aperture
Any experimenter realizing 3D light shaping at some point has to make an implicit or explicit choice about the range of directions of the plane-wave components that superpose in the target volume. This is a choice of the numerical aperture (NA) of the system and it depends on the solid angle spanned by the directions of light rays reaching the target volume. Usually the light-ray shaping. The cone of light angles is determined by the geometry of the setup, in the simplest case the aperture radius of the last lens in the system and its distance from the position where the light is to be shaped (a). The cone angle, α, is also the angle of the k-space-sphere segment (b). The larger the cone angle, the smaller the error ε between the desired and calculated intensity structure (c). The errors in (c) were calculated for the example of a shell-shaped target intensity, shown in the box in the top-right corner. Also shown are the calculated intensity structures corresponding to some of the data points.
directions form a cone of angle α ( figure 5(a)), which in the simplest case is determined by the size of the aperture of the last optical component and the distance between that component and the space where the beam is to be shaped. As each k-space point corresponds to parallel light rays with a specific direction, which is given by the gradient of the phase, such a cone of light-ray directions restricts non-zero k-space values to the section of the k-space hemisphere that lies within an angle α/2 of the k z axis ( figure 5(b)). Figure 5(c) demonstrates the effect of varying the angular size of the k-space-sphere segment on the resulting intensity distributions. It can be seen that a larger cone angle allows the possibility of generating a 3D light intensity that resembles the target intensity more closely. The reason for this is that the larger range of values of k x , k y and k z in such a beam widens the range of structure sizes that can be present in the beam.

Simulation parameters and procedure
Our particular implementation of the 3D GS algorithm, running on a dual 2.5 GHz Pentium Xeon desktop computer, takes several days to converge using a resolution of 256 × 256 × 256. For this reason we initially test the generation of new shapes using a lower resolution of 64 × 64 × 64. In addition to each iteration taking less time at this lower resolution, the algorithm requires fewer iterations to converge.  Figure 6. (a) Schematic of the experimental setup. A collimated HeNe-laser beam illuminates a phase-only spatial light modulator (SLM) [13] before passing through a Fourier lens (L 1 ). The SLM displays a phase pattern that shapes the light beam in a volume around the SLM's Fourier plane. An additional blazed phase grating displayed on the SLM directs the shaped beam into the grating's +1st diffraction order; the other orders-caused by imperfections in the SLM's phase response-are filtered out by a Fourier-plane aperture. The volume in which the light beam is shaped is shown in blue. A CCD, with the help of imaging lenses L 2 and L 3 , records the intensity in a number of planes across this volume. In our experiment f 1 = 600 mm, f 2 = 1350 mm and f 3 = 200 mm. (b) Geometry of the correspondence between position on the SLM and on the k-space hemisphere. Each point light source P in the front focal plane of the Fourier lens L 1 gives rise to a uniform plane wave whose k vector is parallel to the line from P to C, the centre of L 1 .

Experiment
We have created examples of light beams with a shaped 3D intensity structure experimentally, specifically some of the beams shown in figure 2. The experimental setup is outlined in figure 6(a). It creates the shaped beam around the Fourier plane of a phase hologram, which was calculated from the k-space representation of the beam found by the 3D GS algorithm (see section 5.1).
Our hologram patterns were calculated for beams whose k-space representations have a relatively large cone angle of α = 90 • . With our phase hologram, which has a height of approximately 20 mm, this cone angle could be realized by using a Fourier lens with a very short focal length of approximately 10 mm; this would create the shaped beam in a volume of approximate size 0.1 mm × 0.1 mm × 0.1 mm. However, in our experiment we use a Fourier lens with a much longer focal length of f 1 = 600 mm, which corresponds to a significantly smaller cone angle of α ≈ 2 • . Our setup can be seen as a combination of a f = 10 mm Fourier lens that shapes the intensity in a cubic target volume of side length ≈0.1 mm, and two more lenses that image this target volume: a lens with f = −10 mm in the same plane as the 10 mm Fourier lens, and the 600 mm Fourier lens we actually use. Because of the imaging characteristics of this lens pair, the image of the original target volume gets stretched to a size of approximately 6 mm × 6 mm × 500 mm.
In our experiment, an expanded beam from a HeNe laser was reflected off a computercontrolled phase hologram in the form of a phase-only SLM [13]. To deal with imperfections in the SLM, a blazed diffraction grating was added to the phase-hologram pattern. This resulted in the desired beam travelling in the direction of the additional grating's first diffraction order; imperfections in the phase response of the SLM led to additional diffraction orders, which were filtered out by the aperture in the hologram's Fourier plane. We used a lens pair to image a plane at a variable distance z behind the aperture onto a CCD. The transverse magnification was slightly less than one. Intensity cross-sections corresponding to different planes taken across the shaped Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT volume were later combined into volume data, which were visualized with bespoke 3D-viewer software based on the VolPack volume rendering library [14,15].

Calculation of the phase-hologram pattern
The 3D GS algorithm finds phase values for the different points of a monochromatic beam's k-space hemisphere. Points on the k-space sphere correspond to infinite uniform plane waves; the direction of the k vector is the direction of the phase-front normal. Experimentally, a finite uniform plane wave can be created from a point light source in the front focal plane of a lens, whereby the position of the point light source determines the direction of the phase-front normal. Each point in a phase hologram that is placed in the front focal plane of a lens therefore corresponds to a point on the k-space hemisphere and controls its phase. Figure 6(b) shows the geometry of the correspondence between points on the SLM and on the k-space hemisphere. In order to calculate the exact position of the point on the SLM that controls a given point on the k-space hemisphere, the point on the hemisphere should be projected through the centre of the lens into the SLM plane. However, we find that a simpler parallel projection as indicated in figure 1, which gives a phase-hologram pattern that is flipped and distorted, works very well. The distortions are small for points representing plane waves travelling at small angles with respect to the z direction, as is the case in our experiment. We calculate the phase-hologram pattern as the phase of the sum over all k-space elements with the same values of k x and k y . In the limit of an infinitely thin k-space hemisphere, this parallel projection is equivalent to calculating the phase cross-section of the 2D Fourier transform of the shaped beam in the Fourier plane (i.e. at z = 0). Figure 7 shows some of our experimental results. The experimental patterns do not exactly match those calculated by the 3D GS algorithm. We believe this is mostly due to residual astigmatism in the optical system, imperfections in the SLM and imperfect alignment of the intensity data from different planes to 3D volume data. In any case, the results clearly demonstrate that the experiment works.

Conclusions
We have used an algorithm first described in [8] to calculate physically realizable light beams with a shaped intensity distribution, and we have created examples of such light beams in an experiment using a phase hologram calculated from the beam's Fourier-space representation. We are currently considering the use of such light beams in optical tweezers and atom optics.