1. Introduction

Since their discovery [1], non-diffracting (ND) light beams have enjoyed great interest, for example in the areas of optical micro-manipulation [2, 3, 4] and medical imaging [5].

It is easy to understand why ND beams do not change shape on propagation [1]. Like all light beams, they can be understood as superpositions of plane waves. On propagation through a distance z, each individual plane-wave component changes phase by k_z∙z (k_z is the wave number in the z direction). Light beams usually change shape on propagation because they consist of plane-wave components with different values of k_z, and which correspondingly change phase relative to one another on propagation. ND beams consist of plane-wave components that all have the same k_z; on propagation, they change phase in exactly the same way, and therefore retain their relative phase, which in turn means that their interference pattern – the beam – does not change.

The wave number k_z is related to the wave numbers in the x and y directions, k_x and k_y, through the equation

For monochromatic beams (to which we restrict our discussion in this paper) a constant value of k_z is associated with transverse wave vector components that satisfy the equation

Equation (2) describes a circle in the k_x-k_y plane. This can be used to create ND light beams as follows (Fig. 1): illuminate a thin ring aperture (an approximation to a circle) in the front focal plane of a lens (we refer to this plane as the Fourier plane), and behind the lens (specifically in the back focal plane, where the beam’s amplitude is the Fourier transform of the amplitude in the Fourier plane) the light will be approximately non-diffracting. (Note that experimentally created light beams are never perfectly non-diffracting; in the setup discussed here this is due to the fact that the intensity in the Fourier plane is a ring of finite width instead of a circle, and also because the aperture of any real Fourier lens is of finite size [6]). This method has been used in various experiments to create different ND beams [7, 8].

An interesting generalisation of ND beams are self-imaging (SI) beams [9]. If a second bright ring is present in the Fourier plane, it creates a second ND beam whose phase changes linearly with z relative to that of the first beam. The beams periodically go in and out of phase – they beat – and the sum of the two beams changes shape periodically. This can be generalised to more than two ND components: all ND components whose longitudinal wave vectors, k_z, satisfy the equation

 

Fig. 1. Experimental generation of ND beams. A laser beam illuminates a thin ring aperture in the front focal plane (plane A) of the Fourier lens. Behind the lens, the beam is non-diffracting. The beam amplitude in the back focal plane of the Fourier lens (plane B) can be easily calculated as the Fourier transform (FT) of the amplitude in plane A. The intensity distributions on the far left and right illustrate the case of a single bright ring (uniform intensity, planar phase) at A, which leads to the simplest Bessel beam behind the lens.

Download Full Size | PDF

suffer the same phase change, Φ₀, during propagation through the distance ∆z. A light beam consisting exclusively of such components will be self-imaging with a period ∆z. SI beams have also been used in optical micro-manipulation [10].

ND and SI beams have been shaped before: for example, the projections-onto-convex-sets algorithm was used to construct ND beams with a cross-section of arbitrarily positioned bright points [11]; diffractive optical elements for the shaping of longitudinally periodic beams were calculated iteratively [12, 13]; SI beams that contain braided vortex lines were constructed mathematically [14]; and the intensity of ND [15] and SI [16, 8] beams was made to take on the shape of an aperture in a double-Fourier-lens setup.

Motivated by our group’s work in optical micro-manipulation, we describe here algorithms for finding physically realisable ND and SI beams whose intensity distributions in one or more planes approximate given arbitrary patterns. Our algorithms are simple modifications of two hologram-calculation algorithms popular in the field of optical micro-manipulation, the Gerchberg-Saxton [17] and Direct-Search [18] algorithms; our modifications make use of the special Fourier-space properties of ND and SI beams. These Fourier-space properties, which ensure that the beam is either ND or SI as well as monochromatic, restrict the shapes the beam can take. We demonstrate the viability of these algorithms with the help of computer simulations and an experiment which demonstrates the creation of light patterns relevant for optical micro-manipulation.

2. Gerchberg-Saxton (GS) algorithm

In this section we show how to use the computationally very efficient Gerchberg-Saxton (GS) algorithm [17] to construct non-diffracting or self-imaging light beams.

Consider the problem of finding a light beam that has the given intensity cross sections I_A(x,y) and I_B(x,y) in the planes A and B, respectively. For simplicity we concentrate here on planes whose light cross sections are related through a single Fourier transform (like in Fig. 1). As the intensity in the two planes is given, the phase distributions in the planes A and B, Φ_A(x,y) and Φ_B(x,y), must be found. The GS algorithm finds good solutions to this problem by the following procedure [19]. Initially chose a light field in plane A with the desired intensity distribution I_A(x,y) and any phase distribution, for example Φ_A(x,y) = 0. The field (phase and intensity) in plane B is completely determined by this; calculate it by taking the Fourier transform of the field at A. The resulting intensity is not usually going to be I_B(x,y); fix this by simply setting the intensity at plane B to I_B(x,y), leaving the phase unaltered. This in turn changes the field in plane A; calculate it by taking the inverse Fourier transform of the field at B. Now the intensity at A is likely to be wrong (i.e. not exactly I_A(x,y)), so fix it by replacing it with I_A(x,y), again without changing the phase. This process of correcting the intensities in alternating planes is repeated until Φ_A(x,y) and Φ_B(x,y) have converged sufficiently. It can be shown that a measure of the difference between the desired and actual intensity distribution decreases monotonically during each iteration [19]. It usually takes a few dozens of iterations for the algorithm to converge sufficiently well.

 

Fig. 2. Simulation of ND (a) and SI (b) light beams, calculated with the GS method, whose intensity cross-section has been shaped to resemble the target intensity shown in (a). In both cases the focal length of the Fourier lens is f = 1m; the radii of the Fourier-space rings are r =10mm (a) and r =6.06mm, 10mm, 12.77mm and 15.05mm (b), resulting in a self-imaging period of 20mm. Shown are intensity cross-sections over a 1mm×1mm area in different transverse planes. The finite width of each ring (the profile is Gaussian with a width of 300μm) limits the distance over which the beams are non-diffracting or diffraction-free to a finite distance; this can be seen in (c).

Download Full Size | PDF

The GS algorithm as described above successively sets the intensity distributions of a beam in two different planes to the desired intensities there. But if the two planes are related through a Fourier transform there is a different way of interpreting this: the GS algorithm then sets the real-space intensity in one plane and the Fourier-space power spectrum of the same beam to the desired ones. It is then straightforward to use this algorithm to shape a ND or SI beam in one plane by setting the beam's real-space representation to the desired shape and ensuring that the beam is ND or SI by setting the Fourier-space power spectrum to be single- or multiple-ring shaped.

Figure 2 demonstrates the simulated propagation of a ND and a SI beam calculated with this method. The beams’ intensity cross-sections (in all planes in the case of the ND beam, in planes z = 0mm and 20mm in the case of the SI beam) are clearly similar to the target pattern. The intensity cross-sections of the SI beam are noticeably more similar to the target pattern than those of the ND beam; this is due to the larger number of rings in k-space representation of the SI beam, which corresponds to more parameters that can be varied to improve the beam’s intensity cross-section, and the fact that in the particular example shown in Fig. 2 the largest k-space rings of the SI beam are larger than that of the ND beam, which means that higher spatial frequencies and correspondingly more small-scale details are present in the SI beam. The intensity cross-sections of the ND beam are shown over a relatively large propagation range to demonstrate the fact that the beam is ND over a finite range only. On a 600MHz Apple G3 laptop, our program – written in LabVIEW and using a grid of 512×512 pixels – performs one iteration in approximately 5 seconds and takes about 30 iterations to reach good solutions.

A possible extension of this method is the use of the multi-plane GS algorithm [20] or even a full 3D GS algorithm [21, 22], which would allow SI beams to be shaped in multiple planes or even over the entire volume of one period.

3. Direct-Search (DS) algorithm

DS algorithms [18] provide perhaps the simplest way to design light beams. The basic DS algorithm can be summarised as iteration of the following procedure. Make a random change to the light beam in one plane. Model the effect this has on the beam and evaluate its “merit”, which has to be defined such that it indicates to which degree the new light beam has the desired properties. For example, the merit function could be defined such that it is greatest if the new light beam’s intensity profile resembles most closely a pre-defined target profile. If the change has resulted in a better beam keep the change, otherwise discard it. As the merit can only increase, this algorithm has to converge. However, it does not always converge to the global merit maximum, but to a local merit maximum. Fast computers have enabled the useful application of this algorithm.

We apply this algorithm by randomly altering the phase and/or intensity at a fixed number of “source points” on one or more centered circles in the Fourier plane. The remainder of the beam in that plane is dark; as before, this ensures that the resulting beam is non-diffracting (one circle) or self-imaging (several circles). The beam is evaluated at a fixed number of “merit points”, which are points at which the intensity is of particular interest. In optical micro-manipulation, for example, the aim of light shaping is often to achieve high intensity at the positions where particles are to be trapped (and low intensity everywhere else). In this case the trapping positions could be the merit points, and the merit function could be defined such that it is increases if the intensity in the merit points increases.

Here we define the merit of a beam in terms of the intensities, I_i, at the merit points as a fraction of the overall power in the beam, P, as

where n_i is the number of merit points in the same plane and ε is a constant much smaller than the other terms that prevents individual terms from becoming -∞. We choose the form (4) because it prefers an equal distribution of the intensity between all the merit points to simply putting all the intensity into a few or even one merit point (an infinitesimal increase in the intensity leads to a merit increase proportional to the increase divided by the square of the intensity). Each term (n_iI_i)/P describes the fraction of the power ideally in a point (P/n_i) and the intensity there. In this way we preferentially brighten up the darkest merit points in each plane. Simple modifications of the merit function could lead to arbitrary relative merit-point intensities, or even a preference for darkness (no intensity) at some points. Other merit functions could be used to shape the beam in whole areas or volumes, for example by densely (of the order of the wavelength apart) covering the areas or volumes with merit points.

The intensity at the merit points due to the field at the source points can be calculated in a number of ways. We use here the fact that each source point (like every point in the Fourier plane) gives rise to a uniform plane wave behind the Fourier lens, whose complex amplitude is proportional to the complex amplitude at the source point. The overall complex amplitude at a merit point is simply the sum of the complex amplitudes due to all the source points. As only one source point is modified during each iteration, we calculate only the change of the merit-point fields due to the modification of the point source (for each merit point, we calculate the field due to the altered source point alone twice, respectively using the old and new source amplitude; the field change is then the difference between new and old field value). Running on the same 600MHz Apple G3 laptop, our program (again written in LabVIEW), when using 10 merit points, performs about 270 iterations per second, whereby the time per iteration is roughly proportional to the number of merit points. It takes typically of the order of 50 “visits” to each source pixel to arrive at a good hologram, and as we typically use around 500 source pixels this takes roughly 1.5 minutes.

 

Fig. 3. Simulated shaped ND and SI beams calculated with our DS method. (a), (b): intensity cross-sections of ND beams shaped into the constellation “Orion” (a diagram of the star position is inset). (c), (d): intensity cross-sections over two self-imaging periods (∆z = 40mm) of SI beams shaped into stretched unit cells of face-centred cubic (fcc, (c)) and body-centred cubic (bcc, (d)) “crystals” of bright spots; the merit points were located in the planes z = 0 and z = 10mm. In (a) both phase and intensity were modulated, in (b) and (d) only the phase (uniform intensity), in (c) only the intensity (uniform phase). In (a) and (b) r = 25mm (1 circle in the Fourier plane), in (c) and (d) r = 6.06mm, 10mm, 12.77mm and 15.05mm (4 circles); each circle contained 512 source points. In all cases f = 1m and λ = 633nm. All intensity cross-sections represent an area of approximately 0.5mm × 0.5mm. The additional multimedia material contains movies showing a computer rendition of the three-dimensional intensity distributions over slightly more than one self-imaging period of the fcc beam (308K) and the bcc beam (348K) as seen from different viewing angles. [Media 1] [Media 2]

Download Full Size | PDF

Figure 3 shows some examples of intensity distributions calculated with our DS algorithm. In two examples the algorithm modulated only the phase of source points, in one example it modulated only the intensity, and in another example both. These examples demonstrate the creation of arbitrary non-diffracting point patterns, in our case in the shape of the stars in the constellation Orion, and “light crystals”: periodic light distributions in the shape of a series of simple crystallographic unit cells. The detailed intensity is not exactly the desired point pattern, but this can be improved with larger and/or more rings in Fourier space.

4. Experiment

We have performed an experiment in which we demonstrate the holographic generation of a SI beam that has been shaped in more than one plane (shaped ND beams and SI beams shaped in one plane have been created before experimentally [15, 8]). Specifically, we have used an intensity-modulating spatial light modulator to convert a Gaussian laser beam into the bcc light crystal calculated in the previous section (Fig. 3(d)).

To convert one laser beam into another both phase and intensity of the original beam has to be altered (unless either the intensity or phase structure of the original beam happens to be that of the desired beam, which is not the case here). Our experiment can be understood by thinking of the original beam as the sum of the desired beam (travelling at an angle α with respect to the original beam, which is taken to travel along the z axis) and a rest; our setup simply subtracts this rest from the original beam. Specifically, we use essentially a uniform plane wave (in fact a widened and collimated Gaussian beam from a commercial HeNe laser) as our original beam, and pass it through an intensity hologram (in the form of an SLM), which absorbs one part of the rest and transmits the sum of three beams: u(x,y)exp(iksin(α)x), the desired beam (u(x,y)), travelling at an angle α in the xz plane with respect to the original beam; u ^*(x,y)exp(-iksin(α)x), a beam (specifically the complex conjugate of the desired beam) travelling at an angle -α with respect to the original beam in the xz plane; and u ₀ = min_x,y(u(x,y)exp(iksin(α)x) + u ^*(x,y)exp(-iksin(α)x)), a uniform plane wave travelling in the z direction. This particular sum of beams is chosen because it has planar phase fronts, and its generation from a beam with planar phase fronts (like our collimated Gaussian) therefore only requires intensity modulation into the form

The exp(±iksin(α)x) terms give the intensity modulation some characteristics of an intensity grating: the three transmitted beams at angles +α, 0 and -α with respect to the original beam can be seen as the +1st, 0th and -1st diffraction orders in the Fourier plane, respectively. The 0th and -1st orders are subsequently filtered out (Fig. 4(a)).

As the desired beam already has to be present in the illuminating beam, the holographic conversion works most efficiently in a plane where the intensity of the original and desired beams (the latter at an angle) overlap most. To achieve a good intensity overlap with the uniform illumination beam, we avoid generating the desired beam in a plane where its intensity distribution is concentrated in only a few bright rings (like in the image of the hologram plane) or spots, like for example in the planes z = 0 and z = 10mm in Fig. 3(d), and instead generate it in a plane half-way between these two planes, that is by choosing u(x,y) to be the amplitude distribution of the desired beam at z = 5mm. The corresponding intensity-hologram pattern, calculated with equation (5), is also shown in Fig. 4(a).

 

Fig. 4. Experimental setup for holographic generation of shaped self-imaging beams (a) and observed intensity distributions over an area of 3.2mm×3.2mm ((b); the corresponding modelled intensity distributions are those shown in Fig. 3(d)). Insets in (a) show the intensity hologram (corresponding to an area of 8.5mm×8.5mm) and the resulting diffraction orders in the Fourier plane, A; z ₀ = 100mm is the z coordinate of plane B, into which the hologram plane, H, is imaged (after Fourier filtering). We use f ₁ = 119mm and f ₂ = 550mm. L ₂ is the Fourier lens from Fig. 1.

Download Full Size | PDF

Figure 4(b) shows the intensity distributions recorded over one self-imaging period. These experimental intensity distributions are not perfect reproductions of the modelled ones shown in Fig. 3(d); we believe this is mainly due to the difficulty of our specific SLM [23] to display high spatial frequencies. However, the intensity distributions clearly show the desired stretched body-centred cubic (bcc) structure.

5. Conclusions

We have demonstrated that the special Fourier-space properties of ND and SI beams can be utilised in well-known algorithms to design such beams. Specifically we demonstrate, by computer-modelling and experimentally, the generation of SI beams shaped in more than one plane. Because they allow the creation of optical lattices with a wide range of symmetries, we believe that such beams offer interesting new possibilities in research areas such as study of cold atomic gases in optical lattices [24].

Future work could include a more detailed analysis. Specifically, the two algorithms should be applied to the same examples and compared in terms of the merit of their output, which could for example be measured in terms of merit function such as the one given in equation (4) or simply the combined power in all the target points, but also in terms of their speed and suitability to various design problems. The dependence on the number of rings could also be quantified.