## Abstract

Gaussian-apodized Bessel beams can be used to create a Bessel-like axial line focus at a distance from the focusing lens. For many applications it is desirable to create an axial intensity profile that is uniform along the Bessel zone. In this article, we show that this can be accomplished through phase-only shaping of the wavefront in the far field where the beam has an annular ring structure with a Gaussian cross section. We use a one-dimensional transform to map the radial input field to the axial Bessel field and then optimized the axial intensity with a Gerchberg-Saxton algorithm. By separating out the quadratic portion of the shaping phase the algorithm converges more rapidly.

© 2013 OSA

## 1. Introduction

Bessel beams have been explored for a number of years as a means to localize the intensity of a laser beam over a distance much longer than the Rayleigh diffraction range corresponding to the focal spot size. Since an ideal conical wave with a constant approach angle to the optical axis will have the same spot size independent of propagation distance, these beams are often called diffraction-free. Moreover, in principle, linear superpositions of Bessel beams can be used to engineer beams with desired profiles [1]. However, ideal Bessel beams require an infinite amount of energy and thus in practice Bessel beams are realized through truncation of the ideal wave. These apodized Bessel-beams accurately approximate the ideal Bessel beam in a “Bessel-zone” of finite width [2] and transition to a ring-beam in the far field [3, 4].

There are numerous applications of Bessel beams in a wide range of fields. On low energy systems, Bessel beams have been used for large depth of focus microscopic imaging systems [5, 6], optical tweezers [7], and optical coherence tomography [8]. They also have applications in laser micro-machining [9], laser electron acceleration [10], and filament formation [11]. In addition, to the best of our knowledge, the Bessel beam is the most efficient means of producing a high-intensity line focus, much more efficient than a cylindrical lens focus. The Bessel-Gauss beams in particular are useful for projecting a line focus to a position at some distance from the optic, especially for intense beams.

The two common methods for creating apodized Bessel beams are diffractive optical elements and
axicons. The advantage of diffractive elements is that there is more flexibility in engineering the
beam properties, but they are lossy and exhibit strong chromatic dispersion. Refractive and
reflective axicons are more energy efficient but are difficult to fabricate and the Bessel-zone of
such beams is localized near the tip of the axicon. In particular, with uniform input, the axial
intensity profile from an axicon grows linearly with distance near the tip; with Gaussian beam
input, the output axial profile has a nominal profile of *z*exp
(−*z*^{2}/*w*^{2}). As will be shown below,
the axial profile of a Bessel-Gauss focus has a Gaussian shape. For many applications it is
advantageous to create an axial intensity profile that is nearly uniform in the Bessel-zone. The
nonlinear material interactions necessary for micro-machining or plasma filament formation in air
can take place within a sharp threshold and a uniform intensity profile will lead to well-controlled
uniform plasma and the associated material modification. A number of methods have been proposed to
engineer the radial phase profile of the optic to produce a uniform line focus. One of these is the
logarithmic axicon which has been used to generate super-Gaussian profiles on the optical axis
[12, 13]. Improvements in aspheric optical element production techniques have allowed the
fabrication of refractive versions of the logarithmic axicon for constant intensity profile
[13] and also for a profile that produces
constant axial intensity in the presence of linear absorption [14]. A two-step method to first convert a Gaussian beam to a ring
beam with a controlled amplitude profile, then to redirect that beam to a flat top axial line focus
was proposed by Honkanen and Turunen [15].
Fresnel refractive axicons have also been attempted [16]. A second technique to control the axial profile is with holographic diffractive
elements. The first demonstration of this type was by Sochacki et al [17]; the holographic technique has recently been improved by using a
liquid-crystal spatial light modulator (SLM) to create the desired axial profiles [18]. While the holographic method can produce high-fidelity
profiles, they exhibit low efficiency (∼5–7%), since the fidelity in the
target region is achieved by allowing power to scatter away from the area of interest [19]. It is also important to note that in the logarithmic
axicon design the profile of the beam on the optical axis has the same functional form as the radial
cross-section at the front of the axicon. That is, ring-beams with a Gaussian or super-Gaussian
cross section passed through the axicon will form Gaussian or super-Gaussian profiles on the optical
axis [12].

In this paper we present an efficient method for controlling the axial intensity profile of apodized Bessel-beams through a single step of phase-only shaping. The method we present in this paper is more general than the logarithmic axicon approach in that it allows different initial beam radial profiles to be shaped into a desired on-axis profile. Specifically, we consider starting with a Gaussian ring-beam and adjusting its wavefront to produce a desired axial intensity profile. Such a ring beam can be formed, for example, by passing a Gaussian beam through a diffraction grating with concentric grooves [20]. While both phase and amplitude shaping will perform the task, it is either lossy or requires two separate shaping steps. Our phase shaping scheme is inspired by work that has been performed to shape the focal spot of conventional Gaussian beams. An example is shaping the wavefront of a Gaussian beam to create a super-Gaussian beam profile [21] at the focal plane of a lens. The technique we use in this paper to determine the optimum wavefront is the iterative Gerchberg-Saxton (GS) algorithm [22]. The GS algorithm was originally introduced for phase retrieval in electron microscopes [22] and has since found applications in astronomy [23], radial beam shaping [24], pulse shaping [25] and image reconstruction [26]. While it is well-known that the focal plane of a lens is the Fourier transform of the input field, it is less appreciated that there is a Fourier-like mapping of the input radial profile to the axial profile. This allows us to reduce the two dimensional problem to one. In this work we employ a modified version of the GS algorithm with rapid FFT iteration between two domains to optimize the axial profile to a super-Gaussian. We find that the phase shaping can be simplified by controlling the initial divergence of the input beam before the ring Gaussian is formed.

## 2. Bessel-Gauss beams

As mentioned in the Introduction, all physical Bessel beams must be of finite extent. One
attractive representation of a finite conical wave is the Gaussian-apodized Bessel beam, also
commonly called the Bessel-Gauss (BG) beam [3]. By starting with an expression for a Gaussian beam tilted with an angle
*γ*_{0} to the *z*-axis it has been shown within the
paraxial approximation that these beams can be expressed analytically by:

*E*

_{0}is the peak field at the origin (

*z*= 0),

*k*is the wavenumber in free space,

*β*= (

*ω/c*)sin

*γ*

_{0}is the projection of the wavenumber on the optical axis,

*w*

_{0}is the characteristic width of the input Gaussian and we have suppressed the exp(

*ikz*−

*iω*

_{0}

*t*) carrier wave. Departing slightly from the notation used in [3], we have expressed the Gaussian beamlet profile in terms of the familiar Gouy phase Φ(

*z*) and complex beam parameter

*q*(

*z*), which in turn is defined as below in terms of the

*z*–dependent beam radius

*w*(

*z*) and radius of curvature

*R*(

*z*):

*r*= 0 in Eq. (1)) intensity profile takes the form

*w*is constant and 2

*w*/ sin

*γ*

_{0}is the axial projection of the spot diameter. The opposite limit is where there is considerable diffraction of the Gaussian beamlet. In this case, if there is to be a ring beam in the far field,

*γ*

_{0}must be at least twice the divergence angle of the beamlet,

*w*

_{0}/

*z*. In this case, the axial intensity profile can be well-approximated by Eq. (3) if we set

_{R}*w*(

*z*) =

*w*

_{0}.

Bagini et al extended the analytic formulation of the BG beams to create “generalized” Bessel-Gauss beams [4]. In particular, it was shown that the profile given by Eq. (1) can be created remotely by propagating a ring-beam of the form

*I*

_{0}is the zero order modified Bessel function of the first kind and

*a*is the the radius of the ring-beam. The Bessel-zone for these beams is centered at the point

*z*=

_{d}*ak/β*and the width of the Bessel-zone can be accurately approximated by

*z*≈

_{w}*kw*

_{0}/

*β*. Moreoever, if

*a/w*

_{0}≫ 1 the initial intensity is well approximated by

## 3. Phase-only shaping of axial intensity profile

#### 3.1. Mapping *E*(*r*, 0) to *E*(0,
*z*)

As seen above, a Gaussian ring-beam focused with a lens or axicon will produce a BG beam with a
Gaussian axial intensity profile. To control the shape of the axial intensity profile we will make
modifications to the radial input field. Rather than computing the full Fresnel transform, the input
field *E*(*r*, 0) can be transformed directly to the axial field
*E*(0, *z*) [27]. (Note that from this point forward we will treat *z* = 0
as the input plane). Under the assumption of azimuthal symmetry the Fresnel propagation integral
reduces to

*s*=

*r*

^{2}then

*=*

_{z}*k*/2

*z*and

*G*(Ω

*) =*

_{z}*F*(

*k*/(2Ω

*)). Equation 7 is in the form of a Fourier transform with Ω*

_{z}*taking on the role of the frequency variable which allows for rapid numerical evaluation using the fast Fourier transform (FFT).*

_{z}While the mapping from the radial input field to the axial output field described above does not
make any assumptions about the propagation distance *z* except that the Fresnel
transform may be used, we are especially interested in the regime where the distance between the
Bessel-zone is large compared to the (effective) Rayleigh range of the envelope on the Bessel
function in the center of the Bessel-zone. Therefore, we will consider an input ring-beam of modulus
*E*_{0}*f*(*r*) with a modified radial phase
factor *ϕ*(*r*^{2}) focused by a lens of focal length
*z _{d}*:

*a*and characteristic width

*w*

_{0}and by Eq. (5) adequately approximates Eq. (4) when

*a*≫

*w*

_{in}.

Note that the lens phase can also be written in Ω* _{z}*-space as
exp
[−

*i*Ω

_{zd}

*s*], where Ω

*=*

_{zd}*k*/2

*z*. Therefore we can separate the lens phase from the input function,

_{d}*g*(

*s*) and the transform of the ring-beam with the shaper phase

*ϕ*(

*s*) will be a function of the shifted conjugate variable

*′*= 0. Finally, we stress that even though a lens may be used, and there is a Fourier-transform-like relationship between the starting and target domains, the mapping we are using is correct to within the Fresnel, rather than the Fraunhofer limit.

_{z}#### 3.2. Optimization of the input radial phase with the Gerchberg-Saxton algorithm

To construct a phase *ϕ*(*r*^{2}) that will
generate a desired on-axis intensity
|*H _{T}*(Ω′

*)|*

_{z}^{2}we consider the following variational problem

*H*(Ω

*′*) = ℱ[

_{z}*E*

_{0}(

*h*(

*s*)exp(

*iϕ*(

*s*)))] and ℱ denotes the Fourier transform. In phase retrieval, the minimum of (13) is zero since it is assumed that there exists

*ϕ*

^{*}such that

*H*and

_{T}*g*(

*s*) =

*h*(

*s*)exp(

*iϕ*

^{*}(

*s*)) are Fourier transform pairs but in our case it is not true in general that the minimum of (13) is indeed zero. However, by the reverse triangle inequality and Plancheral’s theorem it follows that for all

*ϕ*:

Similar to what has been used in pulse shaping applications, for example Rundquist et al
[25], to optimize
*J*[*ϕ*(*s*)] we use a modified
GS algorithm. In this algorithm, the functions *h*(*s*) and
*H _{T}* (Ω

*) are uniformly discretized at the points*

_{z}′*s*and (Ω

_{i}*′*)

_{z}*respectively and following Fienup’s notation [28] the GS routine is*

_{i}Since we are using the GS algorithm as an error reduction algorithm, i.e.
*η _{j}*

_{+1}<

*η*, it follows that the sequence

_{j}*η*will converge to a fixed value [28]. However, this is not the same as convergence of the algorithm and in particular the sequence

_{j}*ϕ*

^{(j)}may not be converging to a phase that is the global minimum of Eq. (13). This lack of convergence manifests itself in the form of stagnation of the algorithm at local minima. If the algorithm gets stuck before getting close to the target function, we use a modification of the GS algorithm called the input-output algorithm [28]. Here a stronger adjustment is applied: |

*H*

^{(j+1)}(Ω

*)| =*

_{s}′*H*(Ω

_{T}*) −*

_{s}′*α*(|

*H*

^{(j)}(Ω

*)| − |*

_{s}′*H*(Ω

_{T}*)|). A typical choice for*

_{s}′*α*is 0.15. We find best convergence if we set

*α*= 0 after the axial function approaches the target. If the optimization converges,

*ϕ*(

*s*) is the required shaper phase.

#### 3.3. Optimization for a super-Gaussian axial profile in the Bessel zone

To create a profile with a uniform intensity we choose the following target field function as a test of our algorithm:

*a/w*

_{0}≫ 1 it follows by Eq. (15) that in order for the the minimum value of

*J*to be as small as possible that which gives us an approximate relationship as to how the various length scales contribute to the mapping between the peak radial and on-axis intensities. In particular, the peak intensity is inversely proportional to the width of the desired target (

*w*), since the available energy is spread out over a volume that is proportional to

_{T}*w*.

_{T}Starting with choices of *a*, *w*_{0}, *k*,
and *z _{d}*, we pick a value for the target width

*w*and run the GS algorithm. Figure 2 shows the results of such an optimization for input conditions

_{T}*z*= 1m,

_{d}*w*= 5mm, and

_{in}*a*= 50mm and a target width of

*w*= 7

_{T}*mm*. For this set of parameters,

*γ*

_{0}≈ 2.9° and the axial width without shaping would be 1.2mm. The optimized axial profile [Fig. 2(a)] is very close to the target super-Gaussian and the error as defined in Eq. (17) is approximately 10

^{−4}. Even more accurate optimization can be obtained with larger

*w*(see below in Fig. 4) The shape of the optimized radial phase for this run is shown in the blue curve of Fig. 2(b). There is a strong quadratic component that will be discussed below.

_{T}To see how the wavefront shaping affects the beam away from the z-axis, we numerically propagated
the input Gaussian ring beam with the calculated optimum phase profile with Fresnel integration.
Figure 3 shows the intensity profile in the Bessel zone. The
Bessel zone radial profile looks as expected for a Gaussian-apodized Bessel beam, except that there
is a taper to the diffraction rings. This taper indicates that the effective approach angle is
z-dependent, *γ*(*z*). For these conditions, the radius to the
first zero varies from 5*μ*m to 6*μ*m across the line
focus. Figure 5(a) shows a cross-sectional lineout of the
intensity profile *rI*(*r*) at a position that is approximately 17mm
after the focal plane. To illustrate that the strong outwardly-decreasing taper is a profile that is
geometrically consistent with a flat top line focus we also plot the function
*rI*(*r*). Resulting from the initial wavefront shaping, the radial
cross-section of the beam evolves towards a profile resembling a flat top sloped away from the
optical axis. Without any shaping, the half width of the axial profile for these conditions would be
1.2mm, substantially smaller than this choice of target width, *w _{T}*
= 7

*mm*. Note that the optimized phase profile is strongly quadratic in the quantity

*r*−

*r*

_{0}(Fig. 2). This indicates a focusing phase, an idea supported by the image in Fig. 5(b). The first-order phase change to the beam is to defocus the ring so that the beamlets are larger where they cross the z-axis. This change in divergence can work with either sign: from the point of view of intense beam propagation, it is desirable to have the ring focus after the Bessel zone. We took advantage of this effect to improve the optimization convergence by estimating the quadratic component of the focusing phase. A straightforward calculation using the Gaussian beam propagation formula for

*w*(

*z*) (Eq. (2)), projected onto the z-axis is used to predict the effective divergence required to spread the beam to the desired size. After re-optimization, the phase correction excursion is considerably smaller as seen in the black dashed curve in Fig. 2.

#### 3.4. Scaling of shaped axial profiles with wavelength and focal conditions

To explore the limits of what shaped widths could be found, we chose to scan through a range of
values for *w _{T}*, using the predicted quadratic phase as the starting phase
for the beam. In Fig. 4, we show a contour map of the
optimized axial profile as a function of the target super-Gaussian width

*w*for the same test conditions as above,

_{T}*z*= 1m,

_{d}*w*= 5mm, and

_{in}*a*= 50mm. The log of the convergence error is shown in the inset. Provided the target width of the superGaussian is sufficiently large, phase-only shaping yields very good results. For small

*w*, there is a overshoot seen at the edge that may be important for some applications. The 3D image on the right illustrates how the decrease in intensity with target width follows the expected inverse dependence.

_{T}Having determined an optimized shaper phase function, how does this function change when the
system parameters are changed? The radial to axial transform in Eq. (6) transforms between the radial variable *s*
= *r*^{2} and its shifted conjugate variable
Ω* _{z}*′ (see Eq.
(12)). For our case we have optimized for an axial profile that is super-Gaussian axial
profile in the variable

*δz*=

*z*−

*z*. Provided that

_{d}*w*/

_{T}*z*≪ 1, Eq. (12) can be written to show that $\delta z=-2{z}_{d}^{2}{{\mathrm{\Omega}}_{z}}^{\prime}/{k}_{0}$. This means that if we have optimized the radial phase profile for a width ΔΩ

_{d}*′, the corresponding length of the Bessel zone in*

_{s}*z*–space can be calculated as ${w}_{T}=2{z}_{d}^{2}\mathrm{\Delta}{\mathrm{\Omega}}_{s}^{\prime}/{k}_{0}$. Thus the same shaper phase profile will result in a similar axial profile if the wavelength or lens focal length are changed. For example, if we change the lens focal length in our example from 1m to 5m, the target width in Fig. 2(b) would change from 20mm to 500mm.

For broadband pulses this scaling leads to the conclusion that the length of the optimized axial
region is proportional to the wavelength. With an axial length of the profile as a function of
wavelength, the bandwidth will vary along the Bessel zone. This effect will be important for pulses
that are sufficiently short that the fractional bandwidth is appreciable. It has been shown that
Bessel beams can be made achromatic provided that the central transverse wavenumber
*k _{T}* =
(

*ω/c*)

*n*(

*ω*)sin

*γ*

_{0}(

*ω*) can be made independent of frequency [29]. This condition is possible if the ring-beam is formed by normal incidence on a diffraction grating, where the diffracted angle follows

*k*= 2

_{T}*π/d*, where

*d*is the groove spacing. Concentric gratings of this form have been developed [20].

## 4. Summary and Discussion

In this paper, we have shown that wavefront shaping of an annular ring-beam with a Gaussian
cross-section can result in a line focus along the optical axis of nearly constant intensity.
Starting from a Fourier transform integral mapping of the radial input profile to the axial profile,
we used a modified Gerchberg-Saxton algorithm to optimize the focused profile using fast
one-dimensional FFT operations. Across the input radial profile, the optimized phase excursions for
our example are approximately 4*π*, well within the phase range of liquid
crystal spatial light modulators. As mentioned in the introduction, the ring beam can be formed with
a concentric diffraction grating, which also allows the radial spectral dispersion to be exploited
to produce an achromatic beam. With the circular diffraction grating, a simple lens or curved mirror
beam expander can be used to control the input quadratic phase of the beamlet before ring formation.
This removes the quadratic component from the shaper, reducing further the phase excursions required
over most of the beam profile. In fact, the residual phase seen in the black dashed curve of Fig. 2(b) is mostly quartic. A well-controlled spherical aberration
term could remove that component as well. It may even be possible to design a multi element lens
with the correct phase profile, as in the lens axicon designed by Burvall et al [30]. At the expense of added complexity, the ring beam can
be formed with another phase plate, in which case a two step shaping process can be used
[15].

We anticipate this approach will find applications ranging from micro-machining to long-range filament formation in air. In the former case, the goal might be to deliver pulses at an intensity that is a prescribed amount over the damage threshold. In the second application, the line focus can be used to efficiently create an ionized column of desired length projected at a long range. We are currently working on the experimental demonstration of the methods described in this paper. The experimental work as well as the transition to applications will be addressed in subsequent publications.

## Acknowledgments

The authors acknowledge funding support under the MURI AFOSR grant FA9550-10-0561 and C.D. acknowledges funding support from AFOSR under the grant FA9550-10-1-0394. The authors wish to thank Ewan Wright for many useful discussions and bringing to our attention reference [4]. C. D. also wishes to thank Daniel Adams for useful input on the convergence of the G. S. algorithms.

## References and links

**1. **M. Zamboni-Rached, “Stationary optical wave fields with arbitrary longitudinal
shape by superposing equal frequency Bessel beams: Frozen Waves,”
Opt. Express **12**(1), 4001 (2004) [CrossRef] [PubMed] .

**2. **T. Graf, D. Christodoulides, M. S. Mills, J. Moloney, S. C. Venkataramani, and E. M. Wright, “Propagation of Gaussian-apodized paraxial beams through
first-order optical systems via complex coordinate transforms and ray transfer
matrices,” J. Opt. Soc. Am. A **29**(9), 1860–1869
(2012) [CrossRef] .

**3. **F. Gori, G. Guattari, and C. Padovani, “Bessel-Gauss beams,” Opt.
Commun. **64**(6), 491–495
(1987) [CrossRef] .

**4. **V. Bagini, F. Frezza, M. Santarsiero, G. Schettini, and G. Spagnolo, “Generalized bessel-gauss beams,”
J. Mod. Opt. **43**(6), 1155–1166
(1996).

**5. **P. Dufour, M. Piché, Y. De Koninck, and N. McCarthy, “Two- photon excitation fluorescence microscopy with a high
depth of field using an axicon,” Appl. Opt. **45**, 9246–9252 (2006) [CrossRef] [PubMed] .

**6. **T. Planchon, L. Gao, D. E. Milkie, M. W. Davidson, J. A. Galbraith, C. G. Galbraith, and E. Betzig, “Rapid three-dimensional isotropic imaging of living cells
using Bessel beam plane illumination,” Nat. Methods **8**, 417–423 (2011) [CrossRef] [PubMed] .

**7. **J. Arlt, V. Garcés-Chávez, W. Sibbett, and K. Dholakia, “Optical micromanipulation using a Bessel light
beam,” Opt. Commun. **197**4–6),
239–245 (2001) [CrossRef] .

**8. **Z. Ding, H. Ren, Y. Zhao, J. S. Nelson, and Z. Chen, “High-resolution optical coherence tomography over a large
depth range with an axicon lens,” Opt. Lett. **27**(4), 243–245
(2002) [CrossRef] .

**9. **F. Courvoisier, J. Zhang, M. K. Bhuyan, M. Jacquot, and J. M. Dudley, “Applications of femtosecond Bessel beams to laser
ablation,” Appl. Phys. A
(2012).

**10. **B. Hafizi, E. Esarey, and P. Sprangle, “Laser-driven acceleration with Bessel
beams,” Phys. Rev. E **55**(3), 3539–3545
(1997) [CrossRef] .

**11. **P. Polynkin, M. Kolesik, A. Roberts, D. Faccio, P. Di Trapani, and J. Moloney, “Generation of extended plasma channels in air using
femtosecond Bessel beams,” Opt. Express **16**(20), 15733–15740
(2008) [CrossRef] .

**12. **A. T. Friberg, “Stationary-phase analysis of generalized
axicons,” J. Opt. Soc. Am. A **13**(4), 743–750
(1996) [CrossRef] .

**13. **I. Golub, B. Chebbi, D. Shaw, and D. Nowacki, “Characterization of a refractive logarithmic
axicon,” Opt. Lett. **35**(16), 2828–2830
(2010) [CrossRef] [PubMed] .

**14. **I. Golub, T. Mirtchev, J. Nuttall, and D. Shaw, “The taming of absorption: generating a constant intensity
beam in a lossy medium,” Opt. Lett. **37**(13), 2556–2558
(2012) [CrossRef] [PubMed] .

**15. **M. Honkanen and J. Turunen, “Tandem systems for efficient generation of
uniform-axial-intensity Bessel fields,” Opt. Commun. **154**(5), 368375 (1998) [CrossRef] .

**16. **B. Chebbi, I. Golub, and K. Gourley, “Homogenization of on-axis intensity distribution produced
by a Fresnel refractive axicon,” Opt. Commun. **285**(7), 1636–1641
(2012) [CrossRef] .

**17. **J. Sochacki, A. Kołodziejczyk, Z. Jaroszewicz, and S. Bara, “Nonparaxial design of generalized
axicons,” Appl. Opt. **31**(25), 5326–5330
(1992) [CrossRef] [PubMed] .

**18. **T. Čižmár and K. Dholakia, “Tunable Bessel light modes: engineering the axial
propagation,” Opt. Express **17**(18), 15558–15570 (2009) [CrossRef] .

**19. **M. Pasienski and B. DeMarco, “A high-accuracy algorithm for designing arbitrary
holographic atom traps,” Opt. Express ,
**16**(3), 2176–2190
(2008) [CrossRef] [PubMed] .

**20. **L. Niggl, T. Lanzl, and M. Maier, “Properties of Bessel beams generated by periodic gratings
of circular symmetry,” J. Opt. Soc. Am. A **14**(1), 27–33
(1997) [CrossRef] .

**21. **J. A. Hoffnagle and C. M. Jefferson, “Design and performance of a refractive optical system that
converts a Gaussian to a flattop beam,” Appl. Opt. **39**(30), 5488–5499
(2000) [CrossRef] .

**22. **R. W. Gerchberg and W. O. Saxton, “Phase determination from image and diffraction plane
pictures in the electron microscope,” Optik **35**, 237–246
(1972).

**23. **R. A. Gonsalves, “Phase retrieval and diversity in adaptive
optics,” Opt. Eng. **21**, 829–832 (1982) [CrossRef] .

**24. **J. S. Liu and M. R. Taghizadeh, “Iterative algorithm for the design of diffractive phase
elements for laser beam shaping,” Opt. Lett. **27**(16), 1463 (2002) [CrossRef] .

**25. **A. Rundquist, A. Efimov, and D. H. Reitze, “Pulse shaping with the Gerchberg-Saxton
algorithm,” J. Opt. Soc. Am. B **19**(10), 2468–2478
(2002) [CrossRef] .

**26. **J. R. Fienup, “Reconstruction of a complex-valued object from the modulus
of its Fourier transform using a support constraint,” J. Opt. Soc.
Am. A (1987) [CrossRef] .

**27. **A. Lipson, S. G. Lipson, and H. Lipson, *Optical Physics*, 4th ed.
(Cambridge University Press, 2010) [CrossRef] .

**28. **J. R. Fienup, “Phase retrieval algorithms: a
comparison,” Appl. Opt. **21**(15), 2758–2769
(1982) [CrossRef] [PubMed] .

**29. **C. Zapata-Rodríguez, “Ultrafast diffraction of tightly focused waves with
spatiotemporal stabilization,” J. Opt. Soc. Am. B **25**(9), 1449–1457
(2008) [CrossRef] .

**30. **A. Burvall, K. Kolacz, Z. Jaroszewicz, and A. Friberg, “Simple lens axicon,” Appl.
Opt. **43**(25), 4838–4844
(2004) [CrossRef] [PubMed] .